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PREFACE 


The  main  objective  of  this  study  is  to  develop  a  numerical  model 
for  simulating  the  movement  of  well  graded  sediment  through  a  stream 
network.  In  Part  2  of  this  report  the  theoretical  background  of  the 
model  is  described  and  the  governing  equations  are  formulated.  The 
numerical  solution  of  these  equations  is  presented  in  Part  3.,  Hydraulic 
routing  can  be  performed  using  any  acceptable  algorithm  supplied  by  the 
user  because  the  water  movement  is  assumed  to  be  uncoupled  from  the 
sediment  process.  The  model  can  be  used  in  conjunction  with  any 
suitable  sediment  yield  model  to  supply  the  water  and  sediment  runoff 
from  lateral  areas.  The  application  and  results  of  the  model  are 
discussed  in  Part  4.  The  computer  program  for  the  numerical  model  of 
the  East  Fork  River  system  is  described  and  listed  in  the  third 
addendum. 
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square  yards  (sq  yd) 
square  miles  (sq  miles) 
acres  (acre) 
cubic  inches  (cu  in.) 
cubic  feet  (cu  ft) 
cubic  yards  (cu  yd) 
pounds  (lb)  mass 
tons  (ton)  mass 
pound  force  (lbf) 
kilogram  force  (kgf) 
pounds  per  square  foot 
pounds  per  square  inch 
U.S.  gallons  (gal) 
acre-feet  (acre-ft) 


_ To _ 

mi  1 1 imeters (mm) 
meters  (m) 
meters  (m) 
kilometers  (km) 
square  millimeters  (mm2) 
square  meters  (m2) 
square  meters  (m2) 
square  kilometers  (km2) 
hectares  (ha) 
cubic  millimeters  (mm3) 
cubic  meters  (m3) 
cubic  meters  (m3) 
kilograms  (kg) 
kilograms  (kg) 
newtons  (N) 
newtons  (N) 

(psf)  pascals  (Pa) 

(psi)  kilopascals  (kPa) 

liters  (L) 
cubic  meters  (m3) 


Multiply 

_ .bl _ 

25.4 

0.305 

0.914 

1.61 

645 

0.093 

0.836 

2.59 

0.405 

16,400 

0.028 

0.765 

0.453 

907 

4.45 

9.81 

47.9 

6.89 

3.79 


1,233 


1 


INTRODUCTION 


This  report  describes  a  one-dimensional  numerical  model  designed  to 
simulate  sediment  transport  in  natural  channels.  The  physical  processes 
associated  with  the  sediment  movement  are  reproduced  using  a  variety  of 
algorithms.  These  algorithms  incorporate  sets  of  equations  that  operate 
on  input  data  in  a  predetermined  sequence  to  generate  output  data  repro¬ 
ducing  the  actual  physical  process.  A  satisfactory  model  must  include 
the  more  relevant  aspects  of  that  process. 

Sediment  moves  driven  by  hydrodynamic  forces  exerted  by  the  flow  of 
water  which  in  many  instances  is  highly  time  dependent.  The  sediment 
transport  model  must,  therefore,  account  for  unsteadiness  in  sediment 
movement . 

The  dependence  of  sediment  motion  on  flow  conditions  makes  it  also 
dependent  on  the  longitudinal  variations  the  flow  experiences  as  a 
result  of  stream  boundary  irregularities.  These  variations  are 
reflected  in  the  spatial  variability  of  the  sediment  load  distribution. 
Depending  on  particle  size,  some  particles  may  be  carried  primarily  in 
suspension,  while  others  move  entirely  as  bed  load.  In  addition, 
depending  on  flow  conditions,  particles  moving  in  suspension  at  one 
place  may  be  moving  as  bed  load  farther  downstream.  ^  Whether  a 
particular  size  fraction  moves  primarily  as  suspended  load  or  bed  load 
determines  to  what  extent  that  fraction  of  the  sediment  load  will  lag 
behind  the  flood  wave  and  therefore  determines  what  the  magnitude  of 
longitudinal  sorting  will  be.  This  means  that  the  transport  model  must 
reflect  the  dependence  of  the  sediment  load  lag  on  the  material 
properties  of  the  sediment  as  well  as  on  hydraulic  conditions.  Whenever 
the  bed  material  consists  of  a  mixture,  its  transport  involves  the 
motion  of  a  multitude  of  particles  of  diverse  sizes.  Some  particles  may 
deposit  on  the  streambed  while  others  are  scoured  away,  resulting  in  a 
size  composition  of  the  material  in  transport  different  from  that  of  the 
bed.  A  realistic  model  must  account  for  the  interchange  between  the  bed 
surface  material  and  the  moving  sediment  load,  and  should  simulate  the 
residual  transport  capacity  of  the  stream.  The  latter  is  a  measure  of 
the  ability  of  the  flow  to  further  entrain  material  of  a  given  size 
fraction  in  the  presence  of  all  the  fractions  already  in  motion. 


During  the  above  particle  interchange,  the  bed  material  particle 
size  composition  changes  continuously  and,  in  the  process,  the  bed  may 
experience  a  net  amount  of  aggradation  or  degradation.  For  certain  flow 
conditions  the  bed  degradation  in  a  reach  may  be  limited  by  the 
formation  of  an  armoring  layer,  over  which  sediment  may  move  either  in 
suspension  or  as  intermittent  bed-load  waves.  The  armoring  layer  may  be 
destroyed  during  high  flows  and  reformed  at  subsequent  low  flows.  A 
model  must  therefore  be  capable  of  tracking  the  streambed  profile 
evolution  and  the  changes  in  bed  material  size  distribution. 

The  proposed  model  is  designed  to  meet  the  above  criteria.  It  can 
simulate  the  unsteady  transport  of  sediment  mixtures  through  a  network 
of  nonbifurcating  channel  reaches,  and  it  can  be  used  in  tandem  with  a 
suitable  sediment  yield  model,  like  the  one  described  by  Borah  et  al. 
(1981),  which  simulates  the  supply  of  water  and  sediment  runoff  from 
adjacent  upland  areas.  At  present  the  model  is  restricted  to 
consideration  of  noncohesive  bed  materials  only.  It  is  also  restricted 
to  down  channel  streamflow,  and  cannot  consider  the  effect  of  transverse 
currents . 

Parts  2  and  3  of  this  report  provide  a  detailed  description  of  the 
model.  Part  4  discusses  the  validation  of  the  model  on  sets  of 
laboratory  and  field  data.  Coding  details  are  given  in  Addendum  3. 
Some  of  this  material  has  been  presented  in  an  earlier  report  (Borah, 
1979). 
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2  MODEL  FORMULATION 
2.1  EQUATIONS  OF  MOTION 

The  present  model  treats  the  time-dependent  problem  of  one¬ 
dimensional  sediment  routing  in  alluvial  channels.  The  hydraulic 
functions  driving  the  sediment  movement  are  the  local  flow  discharge, 
stage,  and  cross-sectional  flow  area.  Therefore,  an  unsteady  flow 
algorithm  is  required  to  compute  the  time  and  space  distribution  of 
those  flow  parameters  in  the  channel,  given  information  concerning  chan¬ 
nel  geometry,  history  of  inflowing  water  discharge,  and/or  downstream 
stage . 

One-dimensional  hydraulic  routing  algorithms  are  usually  based  on 
the  shallow-water  equations  of  momentum  and  conservation  of  mass  for 
sediment-laden  water.  In  natural  streams  load  concentrations  of  up  to 
50,000  ppm  will  not  change  the  density  of  the  mixture  by  more  than  .3%. 
Thus,  density  variations  may  be  ignored.  Furthermore,  within  that  range 
of  concentrations,  the  surface  waves  propagate  with  a  velocity  that  is 
practically  unaffected  by  the  erodibility  of  the  bed  (Gradowczyk,  1968). 
The  governing  equations  for  water  movement  can  thus  be  written 
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where  A  is  the  flow  cross-sectional  area,  g  is  the  acceleration  of 

gravity,  Q  is  the  flow  discharge,  q  is  the  lateral  water  inflow  per 

w 

unit  length  of  channel,  Sq  is  the  bed  slope,  is  the  friction  slope,  t 
is  t  i.me ,  V  is  the  mean  flow  velocity,  is  the  velocity  component  of 
lateral  inflow  in  the  main  flow  direction,  x  is  the  horizontal  distance 
along  the  channel,  y  is  the  flow  depth,  P  is  the  momentum  coefficient, 
and  p  is  the  water  density  (Fig.  1). 

Two  approximations  to  Eqs.  1  and  2  have  found  wide  application  in 
unsteady  flow  routing.  They  are  the  diffusion  and  kinematic  wave  ap¬ 
proximations  (Ponce,  Li,  and  Simons,  1978).  The  first  assumes  that  the 
inertia  terms  are  negligible,  while  the  second  assumes  that  both  the 
inertia  and  pressure  terms  can  be  neglected.  Both  approximations  have 
neen  shown  to  be  satisfactory  in  a  variety  of  cases.  Solutions  of  the 
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complete  shal  low-water  equations  and  their  v.i ;  n-n.:  ippri.x  i  mat  i  oris  have 
been  extensively  discussed  in  the  literature  (Mi  tier  and  Yevjevich, 
1975)  and  they  will  not  be  considered  in  this  payer  Anv  acceptable 
hydraulic  algorithm  will  suffice  since  the  water  movement  is  assumed  to 
be  uncoupled  from  the  sediment  processes. 

The  third  equation  i  t  motion  is  given  by  the  conservation  of  mass 
equation  for  sediment.  As  shown  in  Addendum  i  ,  this  equation  can  be 
expressed  in  the  form 


9Q 
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9t  3x 
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where  c  is  the  average  volume  concentration.  Q  i  s  the  sediment,  volume 
flux,  A  is  the  bed  poros « *  y ,  W  is  the  active  width  of  the  bed  (i.e., 
that  portion  of  the  fed  width  in  which  erosion  or  deposition  takes 
place),  z  is  the  local  bed  elevation,  r  is  the  net  sediment  volume  flux 
through  the  suspended-bed  load  interface,  and  q  is  the  lateral  inflow 
of  sediment  per  unit  channel  length. 

In  Eq.  3  the  third  term  represents  the  volume  rate  of  sediment 
scour  (or  deposition)  per  unit  length  of  channel  bed.  The  fourth  term 
is  envisioned  as  the  time  rate  of  net  local  exchange  between  the  sus¬ 
pended  and  bed  load  zones .  An  expression  defining  this  term  is  needed 
for  solving  Eq .  3.  For  simplicity,  the  following  parametric  exchange 
equation  is  adapted  (Whitham,  1974) 


3  r 

at 


k  A :  ; 
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.  -It  (R-r)cj 


(4) 


Here  kQ  and  are  constants,  T  is  the  concentration  at  transport  capa¬ 
city,  corresponding  to  which,  k  is  the  net  sediment  exchange  between  the 
suspended  and  bed  zones.  Eq .  4  was  so  constructed  in  order  to  (i) 
incorporate  in  the  solution  some  of  the  nonlinearity  undoubt ly  present 
in  the  exchange  process,  and  (ii)  preserve  the  hyperbolic  character  of 
the  sediment  continuity  equation.  Although  Eqs  .  3  and  4  can  be  solved 
by  successive  iterations  to  obtain  c  and  r,  a  simpler  solution  results 
from  the  following  observation.  Very  near  equilibrium  the  right  hand 
side  of  Eq.  4  vanishes  for  all  practical  purposes.  That  is 


J  .  10 


kQA[(T-c)r-k1(R-r)c]^  0  . 


(5) 


At  the  same  time  T  deviates  slightly  from  its  equilibrium  value.  Thus, 
from  Eq.  5  results 

3r  =  klRT  3c  .  (6) 

at  [T+(k1-l)c]2  at 

Also  this  expression  is  assumed  to  hold  in  slowly  varying  flow  condi¬ 
tions.  Approximating  Qg  by  AVC  over  a  small  time  interval,  and  combin¬ 
ing  Eqs .  3  and  6  gives  the  sediment  continuity  equation  used  in  this 
model 


3c  „  3c  _  qst 

3t  s  8x  A(l+f  ) 

s 

(7a) 

where 

v  =  v 

s  (1+fs}  ’ 

(7b) 

kjRT 

fs  "  AlT+tkj-Dc]2  * 

(7c) 

qst  =  %  '  °-A)W  If 

(7d) 

Eq.  7a  is  a  quasilinear  hyperbolic  equation  governing  the  propagation  of 
the  sediment  concentration  waves.  Its  right  hand  side  represents  the 
lateral  sediment  inflow  contributed  by  runoff,  tributary  channels,  and 
bed  scouring.  The  limited  experience  gained  so  far  with  the  model 

-3 

indicates  that  k^  is  of  order  10  .  Thus  Eq.  7c  has  been  approximated 

by, 


f  £ 
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CEL _  , 

AT[ 1-0. 999c/T] 2 


(7e) 


where  CEL  is  a  user  supplied  parameter.  This  parameter  should  be  ad¬ 
justed  to  match  the  time  of  arrival  of  the  observed  sediment  load  peak 
at  the  channel  outlet.  Eqs.  7b  and  7e  show  that  the  celerity  of  the 
sediment  wave,  V  ,  is  controlled  by  the  local  hydraulic  conditions  and 
sediment  properties  since  capacity  depends  on  these  parameters  as  well. 
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Eq.  7a  can  be  solved  by  means  of  the  method  of  characterist  ics. 
From  this  equation  and  the  total  differential  of  the  sediment 
concentration,  the  following  characteristic  equations  result 


=  V 

dt  s  ’ 

(8a) 

dc  _  ^St 

qst 

•  v  • 

Q  s 

(8b) 

dt  “  A(l+f  ) 
s 

equations  show  that 

Eq.  7a  possesses  only  one 

system  of  forward 

characteristics.  Accordingly,  this  equation  cannot  be  used  in  situa¬ 
tions  where  there  are  flow  reversals.  Integrating  Eqs .  8a  and  8b  with 

the  initial  condition  c  =  c(x  ,t  )  gives 

o  o  o 

x 

c  =  cQ  +  J  (qst/Q)  df ,  (9a) 

X 

o 

x  -1 

t  =  t  +  /  V  (1+1  Jdt  .  (9b) 

o  s 

X 

0 

These  integrals  are  used  to  track  the  evolution  of  the  sediment  waves 
across  the  characteristic  plane.  The  concentrations  existing  on  all  the 
characteristics  at  the  time  of  their  arrival  to  the  downstream  boundary 
define  the  outflow  sedimentgraph . 

It  is  interesting  to  note  that  the  preceeding  equations  reflect  the 
expected  behavior.  For  instance,  it  is  a  recognized  fact  that  the 
streamwise  velocity  of  sediment  particles  always  lags  behind  the  velo¬ 
city  of  the  surrounding  fluid  (Francis,  1971).  This  velocity  lag  ranges 
from  very  small  values  for  silt  particles  in  suspension,  to  quite  large 
differences  in  the  case  of  coarse  sands  and  gravels.  Thus,  the  celerity 
of  a  sed iment- 1 oad  wave  will  be  smaller,  in  general,  than  the  celerity 
ol  the  carrying  flow  wave.  This  trend  is  also  predicted  by  Eqs.  7b  and 
7e,  for  they  indicate  that  waves  of  coarse  material  will  travel  slower 
than  waves  of  finer  material  given  th.it  the  carrying  capacity  of  a 
stream  increases  as  the  sediment  size  decreases. 

This  celerity  lag  is  strictly  a  function  of  local  flow  and  sediment 
properties,  and  it  should  not  he  contused  with  the  differences  sometimes 
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observed  between  the  time  of  arrival  of  the  flow  and  sediment  load 
peaks.  In  fact,  the  sediment  load  may  peak  ahead,  in  phase,  or  behind 
the  flow  peak  depending  on  antecedent  conditions,  intensity  of  the 
event,  sediment  source  location,  season  of  the  year,  etc.  (Guy,  1970). 
To  illustrate  this  point  consider  a  channel  reach  having  a  deposit  of 
very  fine  sediment  on  the  upstream  portion  of  the  reach,  and  receiving 
an  inflow  as  shown  in  Fig.  2a.  The  sediment  will  be  entrained  as 
soon  as  reaches  a  sufficient  intensity,  and  will  move  along  charac¬ 
teristics  parallel  to  the  flow  characteristics  (i.e.,  Vg  =  V,  Eq.  8a). 
As  the  characteristics  reach  the  downstream  boundary  the  water  and 
sediment  outflows  begin  to  rise,  generally  at  different  rates.  The 
concentration  increase  rate  is  controlled  by  the  ratio  of  sediment 
entrainment  to  flow  discharge  (Eq.  8b).  If  the  supply  of  loose  sediment 
is  finite,  the  sediment  outflow  will  peak  at  some  time  t^,  while  the 
water  outflow  will  continue  to  increase  and  peak  at  a  later  time  t^. 
The  lag  between  these  two  peaks  will  obviously  depend  on  the  magnitude 
of  the  sediment  supply,  inflow  rates,  and  distance  of  sediment  source  to 
channel  outlet  (x  =  ,  Fig.  2a).  For  simplicity,  this  example  has 
ignored  backwater  effects  to  restrict  the  flow  movement  to  forward 
characteristics.  Next,  consider  the  same  channel  reach  but  with  the 
source  located  farther  upstream  and  a  constant  base  flow  0^  (Fig.  2b). 
The  flows  and  sediment  loads  reaching  the  channel  inlet  (x  =  x^)  are 
generated  by  characteristics  emanating  from  the  upstream  reach,  xq  ^  x  ^ 
Xj.  Lateral  runoff,  q^,  is  assumed  to  begin  entering  the  channel  at  t  = 
t^.  At  this  point  the  outlet  hydrograph  will  begin  to  rise  and  even¬ 
tually  peak  at  t  =  t,.  .  On  the  other  hand,  the  sediment  characteristics, 
which  enter  the  channel  reach  after  t^ ,  will  reach  the  outlet  later  than 
t^  and  will  finally  peak  at  a  later  time  t^  >  t^.  Thus,  in  this  case 
the  sediment  peak  lags  significantly  behind  the  flow  peak  as  result  of 
changes  in  the  channel  water  inflow  and  sediment  source  location. 

Consider  now  a  channel  with  sediment  moving  mostly  as  bedload.  In 
this  case  the  term  9r/3t  in  Eq.  3  can  be  neglected.  Furthermore, 
Mahmood  (1971))  has  shown  that  ignoring  the  time  derivative  of  the  spa¬ 
tial  concentration  has  little  influence  on  the  simulation  of  bedload 
propagation.  If  in  addition  the  period  of  the  flow  velocity  field  is 
very  large,  the  surface  gravity  waves  can  be  filtered  out  and  only  the 
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2.  Routing  of  water 


sediment 


bedload  propagation  waves  need  be  retained  (Gradowczyk,  1968).  In  such 
an  instance,  Eqs.  1,  2,  and  3  reduce  to  the  equations  governing  the 
"known  discharge"  approximation  frequently  employed  in  modeling  propaga¬ 
tion  of  bed  transients  (Gradowczyk,  1968;  de  Vries,  1971;  Cunge  and 
Perdreau,  1973;  Thomas  and  Prashun,  1977;  Ponce  et  al.,  1979). 

The  foregoing  observations  demonstrate  that  the  equations  presented 
in  this  Section  incorporate  the  ingredients  needed  to  simulate  the 
actual  physical  process.  These  equations  are  complemented  with  a  number 
of  ancillary  algorithms  discussed  in  the  following  sections. 

In  natural  channels  there  is  usually  a  wide  gradation  of  sediment 
sizes.  Although  the  finer  particles  comprise  the  bulk  of  the  load,  the 
channel  evolution  is  controlled  by  the  coarser  particles  (i.e.,  armor¬ 
ing).  Moreover,  different  sizes  are  transported  at  different  rates  as 
pointed  out  earlier.  It  is  thus  important  to  predict  the  movement  of 
the  individual  particle  sizes  encountered  in  the  channel.  The  present 
model  divides  the  size  range  in  a  suitable  number  of  size  fractions,  and 
then  uses  Eqs.  7b,  7d,  7e,  and  9  to  track  the  movement  of  each  fraction. 
2.2  ANCILLARY  ALGORITHMS 

2.2.1  Sediment  Transport  Formulas 

These  formulas  are  used  to  determine  the  potential  carrying 
capacity  or  a  specific  flow.  Different  capacities  can  be  expected  for 
different  particle  sizes,  and  not  all  L.isting  formulas  perform  equally 
well  for  all  sizes.  In  the  present  model  several  formulas  are  used  that 
are  framed  for  easy  use  in  digital  applications,  and  require  only 
information  on  the  hydraulic  parameters  of  the  carrying  flow.  They  are 
the  total  load  formula  of  Yang  (1973)  used  to  estimate  the  transport  of 
very  fine  to  coarse  sands  (0.1-2  mm),  a  duBoys-type  bedload  formula 
(Graf,  1971)  used  with  fine  gravel  particles  (2-4  mm),  and  the  Meyer- 
Peter  and  Muller  bedload  formula  (Meyer-Peter  and  Muller,  1948)  used 
with  particles  in  the  medium  to  very  coarse  gravel  range  (4-65  mm). 
These  formulas  are  presented  in  Addendum  2  in  the  forms  used  in  this 
model.  The  Yang  formula  was  found  to  give  very  reliable  estimates  for 
flows  carrying  sands  in  the  indicated  size  range,  including  flows 
transporting  sediment  mostly  as  bed  load  (Alonso  et  al.,  1980).  The 
duBoys-type  formula  and  the  Meyer-Peter  and  Muller  formula  were  selected 
because  they  were  developed  from  data  in  the  specified  size  ranges. 
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Nevertheless,  the  user  may  replace  these  formulas  by  others  he  may  deem 
more  appropriate  for  a  particular  simulation.  In  particular,  simple 
relationships  between  the  sediment  transport  rate  and  the  flow  condition 
developed  from  in  situ  field  surveys  should  be  given  preference. 

2.2.2  Residual  Capacity.  Composition  of  Material  in  Transport 

Whenever  the  bed  material  consists  of  a  mixture  of  different 
sediment  sizes,  the  particle  size  composition  of  the  sediment  discharge 
usually  differs  from  that  of  the  bed.  Therefore,  the  transport  must  be 
characterized  by  a  load  rate  and  a  size  distribution  analysis.  Although 
a  flow  may  have  the  potential  to  transport  sediment,  its  residual 
capacity  or  ability  to  carry  any  additional  Load  depends  on  the  sediment 
material  already  present  in  the  flow.  Consider,  for  instance,  a  flow 
carrying  a  load  c^  of  uniform  size  d^  and  let  be  the  corresponding 
potential  capacity  of  the  flow.  Then 


Tci  =  Ti  -  = T,  -  Ti  <VTi> 

is  the  residual  capacity  of  the  flow  for  that  size  material.  The  last 
term  in  this  expression  represents  that  portion  of  1  already  consumed 
by  the  material  in  transport.  Similarly,  if  the  later  were  of  a  dif¬ 
ferent  size  d2>  the  residual  capacity  for  the  size  d^  materia.  ,  in  the 
presence  of  the  load  c^ ,  would  ne 

Trl  =  T1  -  T1  'W  ' 

where  c2^2  ma^  env^si°ne^  as  Chat  part  of  T^  depleted  by  the  load 
c^ •  If  both  sizes  were  simultaneously  present  in  the  flow,  then 

Trl  =  T1  '  VW  •  T1  (C2/T2)  ' 

This  expression  can  be  generalized  to  any  size  fraction,  d^  ,  and  for  an 
arbitrary  number,  n,  of  load  fractions  Cj,  c^ ,  ...»  c  ,  as  follows 

Tri  =  Ti  -  WV  -  WV 

-  T.(c./T.)  -  • ' ’  -  T.(cn/Tn),  (10) 

or, 
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The  quantity  within  brackets,  1  -  1  (c./T.),  represents  the  portion  of 

j=i  J  J 

the  potential  capacity  T^  taken  up  by  all  the  size  fractions  in 
transport.  Therefore,  the  quantity  within  brackeLs,  A,  is  the  remaining 
capacity  for  transporting  additional  material  of  size  d^ .  For  a  given 
sediment  load,  A  is  the  same  for  all  fractions.  T.  depends  uniquely  on 
the  local  flow  and  the  properties  of  the  d^  fraction,  while  T  .  depends 
on  all  these  parameters  and  on  the  size  composition  of  the  sediment 
load . 

When  A  >  0,  any  size  fraction  available  for  entrainment  at  the  bed 
surface,  and  for  which  /  0,  will  be  removed  by  the  flow  and  added  to 
the  same  sediment  size  class  already  in  transport.  Thus,  A  >  0  identi¬ 
fies  an  eroding  bed  condition.  Similarly,  when  A  <  0  the  stream  carries 
a  load  in  excess  of  its  potential  capacity  and  will  deposit  the  excess 
sediment  material  on  the  bed.  Therefore,  A  <  0  characterizes  an  ag¬ 
grading  bed  condition.  When  A  =  0  there  is  no  load  change  and  the 
transport  process  remains  in  a  pseudo-equilibrium  condition.  By  contin¬ 
uously  tracking  the  value  of  A,  the  dependence  of  the  individual  resi¬ 
dual  capacities  on  the  composition  of  the  sediment  load  and  its  exchange 
with  the  bed  is  readily  simulated. 

The  size  composition  of  the  total  load  is  continuously  up-dated  by 
adjusting  the  concentration  of  individual  fractions  according  to  the 
composition  of  the  material  removed  from  or  added  to  the  bed.  This 
procedure  is  explained  later  in  the  report. 

2.2.3  Bed  Composition 

In  cohesionless  beds  the  material  available  for  entrainment  is 


essentially  that  exposed  at  the  bed  surface.  As  dunes  and  ripples  move 
slowly  downstream  they  continuously  mix  all  the  sediment  they  contain. 
The  space  occupied  by  those  bed  features  may  thus  be  regarded  as  a 
mixing  zone  below  which  the  bed  material  remains  undisturbed  (Fig.  3a). 
Where  the  sediment  contains  a  mixture  of  different  sizes,  the  slowly 
moving  coarse  material  tends  to  collect  at  the  base  of  the  mixing  zone 
thus  forming  a  lense  of  large  grains.  Under  some  flow  conditions  the 
coarse  material  in  the  bed  may  become  immobile.  In  this  case,  the  flow 
washes  the  finer  particles  out  of  the  mixing  zone  leaving  an  armor  coat 
that  protects  the  underlying  material.  During  the  scouring  of  the  finer 
materials  the  exchange  between  bed  and  flow  takes  place  in  a  thin  layer 
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Fig.  3.  Schematic  representation  of  bed  processes 


at  the  bed  surface,  identified  here  as  the  active  layer.  Several  of 
these  layers  will  be  scoured  away  while  the  mixing  zone  is  degrading. 
When  the  bed  is  armored,  the  last  active  layer  becomes  the  armoring 
coat.  Conversely,  during  the  process  of  bed  aggradation  several  active 
layers  will  be  deposited  on  the  bed  forming  a  new  mixing  zone. 

To  model  the  above  processes  the  mixing  zone  is  pictured  as  a  band 
of  constant  thickness  divided  into  several  layers  (Fig.  3b).  The  layer 
in  contact  with  the  flow  is  always  referred  to  as  the  active  layer.  The 
thickness,  porosity,  and  size  distribution  of  this  layer  can  vary 
throughout  the  simulation,  but  the  layer  is  assumed  to  be  homogeneous 
within  itself  at  any  given  time. 

When  all  the  material  within  the  active  layer  moves,  a  reasonable 
upper  bound  of  its  thickness  is  obtained,  from  volumetric  considera¬ 
tions,  as 


ALT 


100 

P 

n 


(12) 


Here  d^,  P^,  and  A^  are  the  size,  percentage,  and  porosity  of  the 
coarsest  fraction  in  the  sediment  mixture.  For  instance,  in  a  uniform 
bed  material  with  A  =  0.50,  Eq.  12  gives  ALT  =  2d^  which  is  in  agree¬ 
ment  with  the  bed  load  thickness  proposed  by  Einstein  (1950).  However, 
when  the  active  layer  becomes  an  armor  coat,  its  actual  thickness  may  be 
considerably  less  than  that  predicted  by  Eq.  12  because  the  bed  may  be 
armored  by  particles  smaller  than  d  .  In  that  case,  it  is  proposed  to 
compute  the  layer  thickness  from 
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where  d^  is  the  smallest  grain  fraction  the  flow  cannot  transport  (i.e., 
T.=0,  i  =  £,£+1 , . . . ,n) .  Eq.  13  is  also  a  better  measure  of  the  active 
layer  thickness  when  some  of  the  fractions  in  the  bed  layer  cannot  be 
eroded  by  the  flow.  At  low  discharges  only  the  smaller  fractions  will 
be  set  in  motion  and  Eq.  13  will  thus  predict  a  thin  active  layer.  This 
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is  intuitively  correct  since  little  bed  material  will  be  scoured  by  a 
low  flow  during  a  simulation  time  step.  As  the  discharge  increases,  the 
coarser  fractions,  usually  present  in  small  percentages,  will  be  en¬ 
trained  and  Eq.  13  will  predict  a  thicker  layer.  This  behavior  is  in 
agreement  with  the  fact  that  a  greater  depth  of  bed  can  be  sorted 
through  by  a  higher  flow  in  the  same  amount  of  time.  Eq .  13  is  thus 
adopted  as  the  general  expression  for  the  active  layer  thickness.  This 
equation  incorporates  Eq.  12,  in  the  limit,  by  letting  S.  =  n  when  all 
fractions  are  in  motion.  The  sediment  contained  in  the  active  layer  is 
the  only  material  available  for  erosion  during  a  simulation  time  step. 
When  the  bed  is  armored  no  erosion  can  occur  until  the  flow  develops  the 
necessary  T  for  the  smallest  size  fraction  present  in  the  armor  coat. 
When  this  happens  the  armor  coat  becomes  again  an  eroding  active  layer. 
If  deposition  of  a  certain  amount  of  sediment  occurs  during  simulation, 
this  material  is  added  to  the  bed  and  a  new  active  layer  thickness  is 
computed  based  on  the  new  mixture  composition. 

In  order  to  account  for  the  time  evolution  of  the  active  layer 
thickness,  it  is  necessary  to  continuously  track  the  size  composition  of 
this  layer.  This  is  done  by  introducing  the  following  accounting  algo¬ 
rithm.  Consider  a  well  mixed  active  layer  with  three  size  fractions 
C*l<<^2<<^3  eroded  by  a  flow  with  sufficient  transport  capacity  to 
scour  all  three  fractions  (Fig.  4).  The  bed  scouring  is  imagined  to 
begin  with  the  entrainment  of  the  d^-particles  exposed  at  the  bed  sur¬ 
face.  Next,  the  d^-particles  at  the  bed  surface  are  removed,  followed 
by  the  d^-grains  hidden  underneath  the  d^-grains.  Finally,  the  d^-par- 
ticles  are  washed  off  the  bed  surface  followed  by  the  d^  and  d^'parti- 
cles  underneath  them,  and  by  the  d^-grains  uncovered  by  the  removal  of 
the  d^-grains.  Thus,  one  d^-particle  is  removed  for  every  d^-particle 
entriined  by  the  flow,  and  one  d^  and  two  d^-grains  are  associated  with 
the  removal  of  every  dr^-grain.  This  ordering  can  be  readily  extended  to 
»  y  number,  n,  of  size  fractions,  and  il  is  summarized  in  the  following 
"entrainment  frequency"  matrix. 
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Fig.  4.  Schematic  reprcspnt.it  ion  oi  active  layer  composition 
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where  i,  j  =  1,2,3,  ...,  n.  In  this  matrix  the  diagonal  elements  indi¬ 

cate  that  every  size  fraction  at  the  bed  surface  is  scoured  once.  The 
off-diagonal  elements  in  each  row  indicate  the  number  of  times  each 

fraction  d.,  l^i^i,  becomes  available  for  entrainment  once  d.  is  removed 
J  i 

from  the  bed.  Alternatively,  the  elements  of  any  column,  j,  express  the 

number  of  times  the  d^  fraction  is  depleted  due  to  the  entrainment  of  a 
fraction  d^.  Obviously,  the  set  of  entrainment  events  in  Eq.  14  is  one 
of  many  possible  distributions  since  in  actuality  any  number  of  d. -par¬ 
ticles  can  be  associated  with  the  removal  of  a  d^-grain.  In  a  more 
general  stochastic  representation  the  F  matrix  would  be  replaced  by  the 
conditional  probability  matrix  of  the  entrainment  process.  Neverthe¬ 
less,  the  proposed  algorithm  is  adopted  for  its  deterministic  simpli¬ 
city.  It  should  be  noted  that  the  F  matrix  is  time  invariant  and  de¬ 
pends  only  on  the  number  of  fractions  used  in  the  simulation. 

The  amounts  of  eroded  materials  corresponding  to  the  frequencies 

F..  depend  on  the  volumes  of  the  fractions  contained  in  the  active 
ij 

layer.  To  convert  those  frequencies  to  volumes  let  be  the  total 

volume  of  the  d^-size  present  in  the  layer  per  unit  channel  length,  and 

v.  .  the  portion  of  V.  that  becomes  available  for  entrainment  when  the 
ij  J 

d^-fraction,  occupying  the  volume  ,  is  eroded.  Thus, 

v.  .  =  F. .  V.  , 
ij  iJ  i 
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from  which, 
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From  this  expression  one  obtains 
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Eq.  15  can  be  rewritten  as 
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are  the  percentages  of  active-layer  material  in  each  size  class  inter 
va  1  . 


Eq.  17  defines  the'  element:,  of  a  square  matrix  identified  in  here 
as  the  "volume  entrainment  matrix,"  v.  The  order  of  this  matrix  is 
equal  to  the  total  number  of  size  classes  in  the  active  layer.  Its 
diagonal  elements  represent  the  volumes  of  the  individual  fractions  on 
the  bed  surface.  The  off-diagonal  row  elements  contain  the  individual 
fractional  volumes  exposed  bv  the  erosion  of  larger  sizes.  Adding  up 
these  elements  gives  the  vojume  of  potential  erosion  associated  with  the 
removal  of  the  largest  size  in  the  row.  On  the  other  hand,  summing  all 
tne  elements  in  each  column  yields  the  total  volume  of  each  size  class 
present  in  the  active  layer  and  available  for  scour  (Eq.  16).  These 
concepts  are  schematically  illustrated  m  Fig.  5a,  which  depicts  the 
v-matrix  associated  with  a  small  cluster  of  bed  particles  grouped  in 
three  different  size  rlasses.  The  shaded  area  in  Fig.  5b  represents  all 
the  material  that  could  be  entrained  along  with  tile  third  fraction.  The 
shaded  portion  of  Fig.  5c  '(-presents  instead  the  total  volume  of  the 
first  fraction  contained  i  n  the  cluster.  During  bed  degradation,  or 
aggradation,  the  elements  .  I  the-  '.-matrix  are  continuously  adjusted  as 


I 


explained  in  the  next  two  sections.  At  the  end  of  each  simulation  time 
step  the  total  volumes  of  the  individual  fractions  left  in  the  active 
layer  are  introduced  into  Eq.  18  to  calculate  the  size  distribution  of 
the  layer. 

2.2.4  Bed  Erosion.  Armoring 

Whenever  A  >  0  the  bed  is  in  a  degrading  mode.  This  is 
conceptually  modeled  by  depleting  the  elements  of  the  v-matrix,  one  row 
at  a  time,  beginning  from  the  smallest  fraction.  The  diagonal  elements 
are  depleted  first,  for  they  are  exposed  to  the  flow  when  simulation 
starts.  Next,  the  off-diagonal  elements  are  scoured  one  by  one 
beginning  from  the  smallest  size.  When  the  flow  residual  capacity  is 
not  sufficient  to  remove  all  the  >  Junes  in  a  particular  row,  equal 
amounts  (for  simplicity)  are  removed  from  all  the  elements  in  the  row 
until  the  residual  capacity  is  satisfied.  The  eroded  volumes  are  added 
to  the  individual  fractions  already  in  transport,  and  the  value  of  A  is 
recalculated.  This  process  continues  until  either  A  is  no  longer 
positive  definite,  or  the  entire  active  layer  is  worked  through.  For 
example,  the  volumes  (v  ) ,  (^v^,  ^V21  ’  ^V22’  ^V21  ’  etc-^>  ^v33^* 
V3 1  /3 ’  v32,/3;  V33^’  v31/^’  v32 etc-^  would  be  eroded,  in  this  order, 
out  of  the  v-matrix  shown  in  Fig.  5.  However,  the  extent  to  which  these 
volumes  are  actually  entrained  depends  on  the  degree  of  detachabi 1 i ty  of 
the  sediment  particles.  These  concepts  are  used  to  construct  the 
following  algorithm  expressing  the  total  volume  of  sediment  eroded  out 
of  each  size  class  in  the  active  layer  during  a  simulation  time  step: 
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i=l,  2,  n.  In  these  equations  ERO  is  a  user  supplied  erodibility 
parameter.  This  parameter  governs  the  actual  amount  of  bed  material 
available  for  erosion  during  a  simulation  time  step.  ERO  is  calibrated 
by  fitting  the  sediment  yield  volume  to  observed  data.  Everytime  a  new 


e_  is  computed  the  concentration  of  the  corresponding  load  fraction, 
Cj ,  is  upgraded  by  letting 


.  . 

C*  =  C  .  +  — .  ( 20  ) 

J  J  A 


This  concentration  is  then  entered  in  Eq.  10  to  update  A.  When  A  be¬ 
comes  nonpositive  for  any  particular  size  fraction  the  transport  is 
termed  "capacity  limited."  On  the  other  hand,  if  A  remains  positive 
after  depleting  a  fraction  the  situation  is  termed  "supplied  limited." 
In  particular,  if  A  >  0  after  considering  all  fractions  present  in  the 
active  layer,  the  bed  is  said  to  be  armored  and  the  active  layer  becomes 
an  armor  coat.  Fig.  6  summarizes  the  foregoing  simulation  sequence. 

After  the  active  layer  has  been  worked  through  by  the  flowing 
water,  its  volumetric  composition  is  given  by  V^-Ev,  i=2,  2+1,..., m, 

where  2  and  m  are  the  smallest  and  largest  fractions  left  in  the  layer. 
The  actual  thickness  of  this  material  is  then 

m  V .  E* 

AlT*  =  V  *  '  (21) 

1=2  i 

Whenever  ALT"  is  equal,  or  larger,  than  the  thickness  obtained  from  Eq. 

12,  the  later  is  taken  as  the  new  thickness  of  the  active  layer.  In 

some  instances,  however,  ALT*  turns  out  smaller  than  ALT.  In  these 

cases,  a  thickness,  6,  of  undisturbed  material  is  added  to  the  active 

layer  (Fig.  7a).  Because  the  model  does  not  track  the  composition  of 

the  undisturbed  layer,  this  is  always  formed  by  original  bed  material. 

This  material  contains,  in  general,  a  smaller  percentage  of  the  largest 

size  class  than  the  ALT*  layer  (Fig.  7b)  and,  therefore,  6  must  be 

reduced  to  account  for  the  difference.  For  simplicity,  a  simple  linear 

correction  is  assumed  given  by  the  ratio  between  the  percentages,  and 

Pu ,  in  the  eroded  and  undisturbed  layers  (Fig.  7b).  Thus,  the  corrected 
m 

thickness  of  the  new  active  layer  becomes 


J.26 


ALL  SIZES 


J .  2  1 


SIZE  FRACTION  dm 


Fig.  7.  Adjustment  of  active-layer  thickness 
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(22) 


Pa 

ALT  =  ALT*  +  6  . 

c  pu 

m 

The  same  criteria  is  used  when  updating  an  armor  coat  thickness,  except 

that  in  this  case  the  ratio  Pa/Pu  is  replaced  by  P^/P1?. 

mm  r  £  £ 

Finally,  the  local  bed  elevation  of  the  channel  bed  is  updated  by 
subtracting  from  it  the  thickness  of  the  eroded  material,  giving  (Fig. 
7a) 


n  E . 

z  =  z  -  -  1  TV-—-.  •  (23) 

new  old  W  .  ,  ( 1 - A . ) 

i  =  l  i 

2.2.5  Sediment  Deposition 

When  the  model  senses  deposition  (A  <  0),  the  flow  drops  sediment 
on  the  bed  during  the  time  step.  At,  and  the  settled  material  is  added 
to  the  active  layer.  Within  the  present  deterministic  framework, 
deposition  begins  with  the  largest  sediment  fraction  and  continues 
through  the  smaller  fractions  until  either  the  stream  is  no  longer 
overloaded  (A  =  0),  or  all  the  fractions  in  transport  have  been 
depleted.  The  volume  deposited  out  of  any  fraction,  D  .  cannot  exceed 
its  (defect)  residual  capacity,  that  is, 

Maximum  D.  =  T  =  T.  A  .  (24) 

i  ri  i 

However,  whether  this  amount  will  reach  the  bed  during  the  time  step  At 
depends  on  this  being  not  less  than  the  average  time  for  the  sediment 
particles  to  settle  to  the  channel  bed.  Data  by  Jobson  and  Sayre  (1970) 
and  by  Lean  (1971)  indicate  that  this  settling  time  may  be  computed 
using  the  particle  fall  velocities  in  the  quiescent  fluid.  Therefore, 
the  actual  deposition  of  a  size  fraction  during  the  interval  At  is 
calculated  from 
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particles  fall  simultaneously  from  the  water  surface.  The  deposited 
volume  is  subtracted  from  the  material  in  transport  to  yield  the 

size-class  concentration 


From  the  preceding  equations  the  deposition  loop  shown  in  Fig.  8  is 
constructed . 

The  settled  sediment  is  added  to  the  active-layer  fractions,  and 
both  materials  are  assumed  to  mix  thoroughly  yielding  the  volumetric 
composition 


V . ,  when  D .  =  0 , 
1  l 


IV.  +  D. ,  when  D.  *  0, 
i  i  t 

i  =  1,  2,  .  .  .  ,  n.  These  volumes  are  introduced  in  Eq.  17  to  update  the 
bed  composition,  and  this  is  then  used  in  Eq.  12  to  compute  the  new 
active  layer  thickness.  Finally,  the  local  bed  elevation  is  updated  by 
adding  to  it  the  thickness  of  the  settled  material  yielding 
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D,  (Eq.  25) 

r 

C* (Eq.  26) 
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NUMERICAL  SCHEME 


To  run  the  model  in  a  channel  reach  of  length  L,  the  open  domain  [0 
£  x  $.  L]  x  [t  £  0]  is  covered  with  a  rectangular  grid  with  lines  paral¬ 
lel  to  the  x-  and  t-axes.  There  is  a  single  family  of  characteristics 
associated  with  Eq.  7a.  This  leads  to  the  solution  of  an  initial, 
one-point  boundary  value  problem.  Initial  conditions  are  cross-section 
geometry,  bed  elevation,  average  flow  velocities,  flow  depth,  active 
width,  and  bed  size  composition.  Upstream  boundary  conditions  are  time 
histories  of  water  and  sediment  inflow,  and  size  composition  of  sediment 
load.  Let  Ax  and  At  be  the  length  and  time  increments,  respectively, 
separating  the  grid  lines  (Fig.  9).  The  coordinates  of  the  grid  nodes 

are  x  =  mAx  and  t  =  nAt,  m,  n  =  0,  1,2,  ...  The  value  of  any  var- 
m  n 

iable,  say  A,  at  the  node  (x  ,  t  )  is  designated  A  .  Given  the  above 
’  3  m’  n  °  m 

data  specified  along  the  line  t  =  t  ,  the  modi  is  used  to  compute  the 
sediment  discharge,  load  composition,  bed  elevation,  and  active  layer 
composition  at  all  nodes  on  the  next  line  t  =  tn+|-  For  use  in  the 
numerical  scheme,  the  sediment  transport  equations  and  ancillary  algo¬ 
rithms  are  implemented  as  shown  in  the  diagram  of  Fig.  11. 

The  volumes  of  the  sediment  fractions  available  for  erosion  during 

At  are  obtained  from  the  active  layer  thicknesses  on  the  line  t  =  t  as 
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(29) 


These  volumes  are  used  in  Eq.  17  to  compute  the  volume  entrainment 

►  •  n 
matrix  v  . 

m 

The  flow  routing  scheme  supplied  by  the  user  is  used  to  obtain  the 
hydraulic-parameter  values  at  the  nodes  J  and  K  (Fig.  9),  and  their 
averages  are  used  to  compute  an  average  transport  capacity,  T.  The 
concentration  of  each  load  fraction  at  the  end  of  the  characteristic 
path  is  obtained  from  the  following  discrete  form  of  Eq.  9a 
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where  Q  is  the  average  of  QT  and  Q  ,  and  q  is  the  piecewise  uniform 

J  K  St 

lateral  sediment  inflow  over  Ax.  The  amount  of  sediment  entrainment,  or 
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Fig.  11.  Organization  of  numerical  scheme 
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deposition,  and  the  updated  load  concent  ral  i  <>ii  .  c  •.  are  calculated 

using  either  the  erosion  loop  (Fig.  0 j  or  the  deposition  loop  (Fig.  8) 
depending  on  the  sign  of  A.  If  the  character .  st  i  i  of  a  particular 
fraction  does  not  pass  through  the  node  K  (Fig.  9  j ,  the  i on cent  rat  1  on  at 
this  node  is  obtained  by  linear  interpolation  from 
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if  At*  >  At, 


i  t  At*  <  At . 


(  11a.) 


(31b) 


In  these  equations  Ax*  and  At*  are  given  by  discrete  forms  of  Eq.  9b.  A 
similar  correction  is  applied  to  E*  (Eq.  19)  to  determine  the  actual 
amount  of  material  eroded  in  the  interval  At.  The  new  concentrations 
are  used  to  calculate  the  new  volume  load  fractions 


n+1  t  n+1 

i,m+l  Ci,m+1 


(32) 


and  these  are  introduced  into  Eq.  18  to  update  the  load  size  composi¬ 
tion.  The  foregoing  load  updating  sequence  is  summarized  in  Fig.  10. 
Last,  the  new  bed  size  composition  and  elevation  are  computed.  When  a 
pseudo-equilibrium  condition  is  encountered  (A  =  0)  no  bed  updating  is 
necessary.  In  continuous  simulation  the  above  calculations  are  per¬ 
formed  consecutively  for  each  channel  segment  to  update  all  variables 
over  the  length  of  the  channel .  The  process  is  then  repeated  for  each 
time  interval  in  the  simulation  period. 
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MODEL  TESTING 


Tests  were  conducted  to  verify  the  ability  of  the  model  to  simulate 
various  instream  processes.  Published  laboratory  and  field  studies  were 
scanned  for  data  suitable  for  testing  the  algorithms  describing  proces¬ 
ses  of  bed  scour,  armoring,  and  unsteady  transport.  Three  useful  sets 
of  data  were  found.  The  first  set  was  collected  by  Ashida  and  Michiue 
(1971).  They  performed  a  series  of  laboratory  experiments  to  study  the 
effect  of  sediment  gradation  on  channel  armoring.  The  second  data  set 
was  collected  by  Lane  and  Carlson  (1933)  in  the  San  Luis  Valley  canals 
in  south-central  Colorado.  They  made  measurements  on  the  bed  and  bank 
materials  that  formed  the  canals.  The  beds  of  these  canals  have  become 
armored  over  the  years.  The  data  offer  an  excellent  opportunity  for 
testing  the  model  on  a  natural  system.  Finally,  the  model  was  tested 
using  data  collected  by  the  U.  S.  Geological  Survey  on  a  reach  of  the 
East  Fork  River,  Wyoming  (Mahoney  et  al.,  1976).  These  data  offer  the 
opportunity  of  checking  the  model  algorithms  on  an  alluvial  stream  under 
conditions  of  unsteady  flow.  These  tests  are  discussed  below. 

As  is  usual  in  channels  with  material  consisting  of  a  wide  range  of 
grain  sizes,  the  bed  slopes  of  the  channels  used  in  the  above  three 
studies  were  fairly  steep.  In  these  cases  the  kinematic-wave 

approximation  to  the  equations  governing  unsteady  flow  of  water  is 
applicable.  In  this  approximation  the  momentum  equation,  Eq .  1,  becomes 
Q  =  KIN  .  A1/2  P_1/2,  (39) 

where 


KIN  =  C 


1-49  r1/6 
n 


(34) 


is  a  kinematic-wave  parameter.  In  these  relationships  P  and  R  are  the 
wetted  perimeter  and  hydraulic  radius  of  the  channel  cross  section, 
respectively,  C^.  is  the  Chezy  coefficient,  and  n  is  the  Manning 
roughness  factor.  Eqs .  .33  and  34  are  also  valid  when  the  flow  is 
uniform  and  steady.  For  the  purpose  of  simulating  the  East  Fork  River 
data  Eqs.  2,  33,  and  34  were  solved  using  the  kinematic-wave  routing 

scheme  developed  by  Borah  et  al.  (1980). 

The  simulations  required  calibration  of  the  model  parameters  to 
obtain  best-fit  of  model  predictions  to  observations.  One  flow  routing 
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parameter,  KIN,  and  two  sediment  transport  parameters,  CKL  and  ERO ,  were 
available  for  adjustment.  CEL  controls  the  travel  time  of  sediment 
load,  and  ERO  governs  the  total  quantity  of  bed  material  available  for 
entrainment  at  the  bed  active  layer. 

4.1  FLUME  ARMORING  STUDY 

Several  laboratory  experiments  were  carried  out  by  Ashida  and 
Michiue  (1971)  to  study  the  effects  of  bed  armoring  on  channel 
degradation.  Various  sediment  mixtures  were  used  with  different  mean 
diameters  and  geometric  standard  deviations.  These  mixtures  were  placed 
in  a  recirculating  flume  with  a  bed  2.62  ft.  wide  and  6S.S  ft.  long.  In 
each  experiment  a  steady  uniform  flow  was  passed  over  the  bed.  The 
hydraulic  conditions  were  chosen  to  purposely  induce  armoring.  No 
additional  sediment  was  introduced  into  the  stream,  and  the  rate  of 
sediment  collection  in  a  trap  at  the  end  of  the  flume  was  used  as  a 
measure  of  the  total  sediment  discharge. 

Data  from  Ashida  and  Michiue' s  run  No.  2  were  chosen  to  check  the 
performance  of  the  model  in  simulating  the  transport  of  cohesionless 
sediment  mixtures  and  the  formation  of  armoring.  The  hydraulic  condi¬ 
tions  used  in  the  run  are  given  in  Table  1.  In  running  the  model,  the 
flume  was  represented  by  four  reaches  of  equal  length  with  a  point 
inflow  of  water  at  the  upstream  end  of  the  channel.  The  initial 
particle  size  distribution  of  the  bed  material  employed  in  the  run  is 
shown  in  Fig.  12.  In  simulating  sediment  transport,  the  graded  bed 
material  was  represented  by  ten  discrete  particle  size  fractions.  They 
are  listed  in  Table  2.  Fig.  12  shows  the  size  distributions  of  the 
eroded  material  leaving  the  flume  and  of  the  armoring  layer  as  reported 
by  Ashida  and  Michiue,  and  the  best-fit  curves  predicted  by  the  model. 
The  agreement  between  the  simulated  and  the  measured  armored  layer  size 
distributions  is  good.  However,  some  discrepancy  exists  between  the 
curves  for  the  eroded  material.  The  measured  curve  definitely  shows 
larger  sizes  present  in  the  sediment  load.  This  discrepancy  may  be 
attributed  to  inaccuracies  in  the  estimation  of  sediment  entrainment. 
In  the  transport  capacity  formulas  an  average  tractive  force  concept  is 
implied,  which  excludes  the  effect  of  the  large  instantaneous 
fluctuations  associated  with  turbulent  tractive  forces  (Alonso  and 
Coleman,  1981). 
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Table  1  -  Hydraulic  conditions  in  the  flume  and  San  Luis  Valley  Canal 
studies . 


Data 

Source 

Flow 

Rate 

__  (cfs) 

Flow 
Depth 
(ft  ) 

Average 

Ve loci t  y 
(f£s) 

F.nergy 

SI  ope 

Manning' s 
Coefficient 

Ashida  and 
Michiue  (1971) 
Run  2 

1 .06 

0.22 

1.81 

0.00440 

0.017 

Lane  and 

Carlson  (1953) 

128.0 

1  .  77 

4.00 

0.00243 

0.023 
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Table  2  -  Size  distributions  used  for  simulating  the  flume  and  San  Luis 
Valley  Canal  data. 


Flume 

armoring  study 

San  Luis  Valley 

Canal 

Size  Class 

Percent 

Size  class 

Percent 

Interval 

per  class 

Interval 

per  class 

(mm) 

interva 1 

(mm) 

interval 

0.2 

-  0.3 

9.5 

0.000  -  0.149 

4.0 

0.3 

-  0.4 

9.5 

0.  149  -  0.297 

5.7 

0.4 

-  0.6 

11.0 

0.297  -  0.590 

10.1 

0.6 

-  0.8 

6.0 

0  590  -  1.190 

12.0 

0.8 

-  1.0 

4.0 

1.190  -  >.380 

9.2 

1.0 

-  2.0 

12.0 

2.380  -  4.760 

5 . 4 

2.0 

-  4.0 

18.0 

4.760  -  9.525 

13.3 

4.0 

-  6.0 

18.0 

9.525  -  19.05 

17.9 

6.0 

-  8.0 

10.0 

19.05  -  38.10 

13.2 

8.0 

-  10.0 

2.0 

38.10  -  76.20 

9.2 
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fig.  12.  Size  distribution  ci.r/is  oi'i. lined  for  f  J  nine  armoring  test 
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SAN  LUIS  VALLEY  CANAL  TESTS 


The  test  described  in  this  section  is  based  on  data  collected  by 
Lane  and  Carlson  (1953)  on  a  reach  of  one  of  the  canals  located  in  the 
San  Luis  Valley  of  southern  Colorado.  These  canals  were  constructed  in 
the  late  1800's  in  the  alluvial  cone  deposited  by  the  Rio  Grande  River 
on  the  floor  of  the  San  Luis  Valley.  Near  the  appex  of  the  cone  the 
deposits  consist  of  sand,  gravel,  and  cobbles,  the  size  of  the  cobbles 
decreasing  with  the  distance  from  the  appex.  The  canals  were  found  very 
stable,  and  therefore  presented  an  unusually  favorable  condition  for 
studies  of  stable  canals.  Several  reaches  ranging  in  length  from  600  to 
2200  ft.  were  selected  for  study.  Hydraulic  measurements  were  made  on 
the  most  regular  portions  of  the  reaches.  Observations  and  mechanical 
analysis  were  also  made  of  the  bed  and  bank  materials  forming  the  test 
sections.  Samples  of  the  material  in  which  the  canal  was  constructed 
were  obtained  by  excavating  into  the  bank  of  each  section.  Mechanical 
analysis  of  the  material  forming  the  bed  surface  layer  disclosed  that  in 
the  stable  sections  the  finer  material  had  been  removed  from  the  layer 
and  an  armor  coat  of  coarser  material  was  left. 

Data  from  the  test  section  No.  12  was  selected  to  simulate  the 
development  of  the  armor  coat.  The  hydraulic  conditions  observed  in 
this  section  are  presented  in  Table  1.  The  bank  and  bed  material  size 
distributions  are  shown  in  Fig.  13.  For  simulation  purposes,  a  600  ft. 
long  reach  divided  into  four  equal  segments  was  assumed,  with  a  constant 
inflow  of  clear  water.  The  initial  bed-material  size  composition, 
assumed  equal  to  that  of  the  banks,  was  represented  by  the  distribution 
listed  in  Table  2.  The  size  composition  of  the  active  layer  at  2.5, 
12.5,  37.5,  and  87.5  hours  are  shown  in  Fig.  13.  As  time  increases,  the 
layer  is  quite  rapidly  depleted  of  particles  10  mm  in  size  or  smaller, 
and  the  distribution  in  this  size  range  approaches  assymptotically  the 
measured  curve.  However,  the  simulated  distribution  in  the  10  to  75  mm 
r.  i  ze  range  differs  from  the  composition.  This  difference  may  be 
explained  in  part  by  the  use  of  a  constant  rate  of  discharge.  Although 
information  related  to  the  actual  history  of  flows  in  the  canal  is  not 
available,  Lane  and  Carlson  (1953)  mentioned  that  the  canal  had 
sustained  flows  considerably  above  those  observed  dining  their  study. 
Therefore,  the  observed  bed-layer  composition  was  most  probably  shaped 
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Fig.  13.  Evolution  of  the  mzp  distribution  in  to[>  layer  of  San  Luis  Valley 
Cana  1 . 
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4.3  EAST  FORK  RIVER  PROJECT 

A  detailed  program  of  hydraulic  and  sediment  transport  data 
collection  was  undertaken  by  the  U.S.  Geological  Survey  during  1975  on  a 
reach  of  the  East  Fork  River  near  Boulder,  Wyoming.  The  study  reach  is 
approximately  3  miles  in  length  and  has  a  tributary  area  of  about  180 
square  miles.  About  half  of  this  area  lies  within  the  Wind  River 
Mountains.  The  other  half  is  provided  by  the  area  drained  by  the  Muddy 
Creek  tributary.  A  map  of  the  study  reach  showing  the  position  of  the 
principal  gaging  stations  is  given  in  Fig.  14.  The  river  in  the  study 
reach  is  about  60  ft.  wide  and  4  ft.  deep  at  bank  full  stage.  Most  of 
the  water  during  high  flows  comes  from  melting  snow  of  the  mountain 
area.  In  the  East  Fork  River  the  bed  is  gravel  on  the  riffles  and  bars, 
but  coarse  sand  constitutes  the  bulk  of  the  bed  load.  The  bed  material 
in  the  Muddy  Creek  is  sand  but  much  finer  than  the  sand  in  the  East 
Fork.  Stage  recording  gages  were  installed  at  sections  B-l,  B-5M,  and 
B—  1 7 .  Sediment  inflow  to  the  reach  was  measured  at  sections  B-l  and 
B-5M  by  sampling  suspended  sediment  with  a  standard  DH-48  hand  sampler 
and  determining  the  bedload  discharge  with  a  Helley-Smith  bedload 
sampler.  The  bedload  discharge  past  the  section  B- 1 7  was  measured  with 
both  a  bedload  trap  and  a  Helley-Smith  sampler.  Leopold  and  Emmett 
(1976)  and  Mahoney  et  al.  (1976)  describe  in  detail  the  stream,  the 
bedload  trap,  and  the  data  <ollected  during  1975.  The  data  included 
cross-sectional  surveys,  flow  rating  curves,  water  temperature,  sediment 
transport  rates,  and  particle-size  distributions  of  bed  load  material. 
This  data  set  was  also  simulated  by  Bennett  and  Nordin  (1977).  Their 
input  data  reduction  procedures  were  followed  to  some  extent  in  the 
present  test. 

In  this  simulation,  the  study  reach  was  divided  into  the  fourteen 
subreaches  shown  in  Fig.  14.  Surveys  of  cross  sections  were  used  to 
obtain  relationships  between  area,  wetted  perimeter,  and  stage.  Time- 
averaged  active-bed  widths  were  determined  from  the  surveys  by  comparing 
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plots  of  bed  eLevation  at  the  beginning  and  end  ot  the  simulation  per¬ 
iod,  and  observing  which  parts  of  the  bed  moved  vertically.  As  noted  hv 
Bennett  and  Nordin  (1977)  this  prodecure  may  overestimate  the  active 
width  because  not  always  all  parts  of  the  active  width  move  simultan¬ 
eously.  However,  there  may  be  no  better  way,  short  of  continuous  bed 
monitoring.  The  profile  of  six  of  the  cross  sections,  as  measured  on 
June  2,  1975,  are  plotted  i a  Fig.  15.  Also  plotted  are  the  estimated 
active-bed  widths,  and  the  stages  observed  during  June  7.  As  these 
figures  reveal,  overbank  flow  occurred  at  several  sections  along  the 
reach  during  the  high  flow  periods,  including  the  inlet  sections  H- 1  and 
B-5M. 

Flov  input  to  the  system  was  provided  by  point  loads  at  sections 
B-l  and  B-5M.  These  point  loads  were  generated  by  converting  ' !.e  stage 
readings  to  flow  rates  using  the  rating  curves  constructed  from  the 
stages  and  flow  discharges  measured  at  those  sections  The  rating 
curves  reported  by  Mahoney  et  al.  (1976)  are  shown  in  Figs,  lba  and  1 6b . 
The  sum  of  the  inflow  hydrographs  measured  during  the  simulation  period 
is  about  equal  to  the  outflow  hydrograph  observed  at  section  B- 17, 
indicating  that  water  was  carried  by  essentially  translatory  waves. 
This  fact  lends  further  support  to  the  use  of  a  kinematic-wave  routing 
scheme . 

Inflow  of  sediment  at  sections  B-l  and  B-5M  were  obtained  from 
sediment-flow  rating  curves.  These  were  constructed  by  relating  the 
measured  sediment  loads  to  water  discharges  obtained  from  the  above  flow 
rating  curves.  Suspended  sediment  did  not  contribute  signii "cantly  to 
the  total  input  load  and,  therefore,  the  sediment  discharges  were  based 
solely  on  the  bedload  data.  When  plotted  in  log-log  scales  the  data 
exhibited  considerable  scatter  but  displayed  a  linear  trend.  Linear 
interpolation  of  the  logarithmic  values  yielded  the  rating  < urves  shown 
in  Figs.  16c  and  l6d.  The  data  scatter  resulted  in  the  low  correlation 
<  ■•fficients  indicated  in  the  figures.  A  total  ot  ten  different  frac¬ 
tions  with  mean  diameters  0.12,  0 . 2  5 .  0.5.  1.0,  2.0,  4.0,  8.0,  16.0, 
32.0,  and  64.0  mm  were  selected  t  <>  repr<  rut  the  sediment  in  the  bed  and 
in  transport.  The  percentages  >1  inuleii.il  associated  with  each  size 
were  determined  from  the  sampled  i/r  d  ■  ;,t  r  i  but  i  ■••ns  These  percentages 
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were  used  to  convert  the  total  bedloads  to  sediment  transport  rates  for 
each  particle-size  fraction.  Similar  analysis  was  performed  for  the 
data  collected  at  section  B-17.  These  data  are  used  in  comparing 
measured  and  predicted  water  and  sediment  discharge  1 r om  the  study  area. 

A  thirty  minute  time  step  was  used  in  modeling  the  East  Fork  River 
for  a  period  of  twenty-two  days  1 roin  May  29  through  June  19,  197S.  The 
first  step  involved  calibration  of  the  flow  parameter  to  best  fit  the 
observed  outflow  hydrograph.  A  comparison  of  simulated  and  recorded 
flows  at  section  B-17  is  presented  in  Fig.  17.  Agreement  is  satisfac¬ 
tory  at  low  and  intermediate  f I < ws ,  but  the  measured  values  exceed  the 
simulated  values  at  high  flows.  The  reason  for  this  discrepancy  is  the 
occurrence  of  overbank  flows  which  were  not  accounted  for  in  the  con¬ 
struction  of  the  flow  rating  curves  (Figs.  16a  and  1 6  b )  . 

The  next  step  after  obtaining  the  hydraulic  results  was  simulation 
of  the  sediment  transport.  The  simulation  was  performed  a  number  of 
times,  adjusting  each  time  the  sediment  parameters  to  obtain  the  best 
possible  agreement  between  simulated  and  observed  values.  Fig.  18  shows 
predicted  and  recorded  bedload  discharges  at  section  B-17  for  the 
twenty-two  day  simulation  period.  The  dashed  lines  joining  the  data 
points  are  imaginary  lines  used  to  approximate  the  shape  of  the  measured 
sedimentgraph .  The  observed  total  bedload  outflow  between  June  1  and 
June  19  is  1889  tons.  The  predicted  value  is  1942  tons,  less  than  three 
percent  higher.  In  spite  of  this  agreement,  there  is  a  tendency  for  the 
simulated  bedload  rates  to  be  lower  than  the  recorded  values  at  high 
flows,  and  higher  during  the  recessions.  The  differences  are  most 
noticeable  during  the  first  high  flow  period.  A  possible  explanation 
for  these  discrepancies  is  provided  by  field  observation  in  the  East 
Fork  River  recently  reported  by  Meade  et  al.  (1981).  Over  the  years 
they  have  observed  that,  during  the  runoff  season,  some  sections  of  the 
reach  are  scoured  while  other  sections  are  filled,  resulting  in  a 
normni form  storage  of  movable  bed  material  along  the  stream. 
Consequently,  the  relation  of  bedload  rate  to  water  discharge  varies 
significantly  from  one  part  of  the  reach  to  another.  Meade  et  al. 
(1981)  noted  that  i mined i ate  1 v  downstream  of  a  storage  area  the  bedload 
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transport  may  increase  steeply  with  the  initial  increase  in  water 
discharge  and  later,  as  the  stored  movable  material  becomes  depleted  the 
transport  rate  rapidly  decreases,  resulting  in  a  multi-valued  rating 
curve.  They  also  mention  that  this  behavior  has  been  indeed  observed  at 
the  gaging  station  B  —  1 7  (Leopold  and  Emmett,  1976).  This  behavior 
cannot  be  exactly  reproduced  by  the  present  model,  which  assumes  a 
supply  of  transportable  material  uniformly  distributed  along  the 
streambed,  and  single-valued  relations  between  bedload  transport  and 
water  discharge  at  all  sections.  Nevertheless,  the  predicted 
sedimentgraph  displays  the  rising  and  falling  trends  observed  in  the 
recorded  bedload  during  the  twenty-two  day  Lime  span. 

Fig.  19  shows  a  comparison  of  simulated  against  sampled  size  dis¬ 
tribution  of  the  bedload  at  the  beginning,  middle,  and  end  of  the  first 
high  flow  period.  The  model  predicts  fairly  well  the  distribution  of 
the  coarser  materials,  and  the  shift  of  the  d  during  the  event.  The 
recorded  values,  however,  indicate  less  mobility  of  material  in  the 
finer  size  range  than  predicted.  There  could  be  several  reasons  for 
this  disagreement.  For  one  thing  the  kinematic-wave  approximation  used 
in  routing  the  flow  may  be  distorting  the  local  energy  slopes,  resulting 
in  overestimation  of  entrainment  threshold  conditions.  Or,  there  might 
have  been  more  fine  material  being  carried  in  suspension  which  could  not 
be  sampled  as  part  of  the  bedload,  but  which  was  actually  accounted  for 
by  Yang's  total  load  formula  used  in  the  simulation.  A  definite  answer 
requires  further  investigation. 
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5  CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  CONCLUSIONS 

1.  A  one-d imens iona I  numerical  model  has  been  developed  for  simulating 
the  movement  of  well  graded  sediments  through  a  stream  network. 

2.  Hydraulic  routing  is  performed  by  using  any  acceptable  algorithm 
supplied  by  the  user  because  the  water  movement  is  assumed 
uncoupled  from  the  sediment  process.  The  model  can  be  used  in 
conjunction  with  any  suitable  sediment  yield  model  to  supply  the 
water  and  sediment  runoff  from  lateral  areas. 

3.  The  sediment  routing  scheme  is  based  on  the  physical  processes 
governing  the  mechanics  of  sediment  movement  in  alluvial  channels. 
The  model  recognizes  the  effect  of  bed  and  suspended  load 
interaction  on  the  total  load  movement,  can  simulate  bed  armoring, 
changes  in  bed  elevation,  and  longitudinal  sorting  of  eroded 
material . 

4.  The  applicability  of  the  model  is  restricted  to  noncohesive 
materials,  relatively  stable  channel  geometries,  streams  with 
negligible  in  and  out-of-bank  transport,  and  flows  in  which 
transverse  currents  may  be  ignored. 

5.  The  model  gave  satisfactory  results  when  tested  on  laboratory  data 
from  a  flume  armoring  study,  and  field  data  from  the  San  Luis 
Valley  Canal,  Colorado,  and  the  East  Fork  River,  Wyoming.  These 
test  .  tend  to  indicate  that  the  model  adequately  simulates  the 
transport  of  graded  cohesionless  sediments,  including  the  effect  of 
a rmor ing . 

5.2  RECOMMENDATIONS 

i  .  It  is  recommended  that  the  channel  model  he  further  tested  against 

a  variety  of  real  situations  with  special  emphasis  on  the  scour, 
deposition,  and  transport  of  noncohesive  materials. 

2.  it  is  recommended  that  the  model  he  further  developed  and  refined 

to  include  the  following  capabilities:  (a)  improve  the  one¬ 
dimensional  representation  by  separating  flows  in  the  incised 

channel  from  flows  over  flood  plains;  (b)  account  for  in  and  out- 
of-bank  transport,  and  lateral  distribution  of  bed-material 
properties  and  hydraulic  conditions;  (c)  predict  the  variation  of 
lateral  bed  slope  and  lateral  sediment  sorting  around  channel 
bends;  and  (d)  simulate  sediment  retention  by  grasses  and  other 
vegetative  covers. 


J  .  52 


3.  The  channel  model  should  be  tested  using  hypothetical  situations  to 
confirm  that  the  model  does  respond  in  a  realistic  manner.  for 
instance,  these  tests  may  include  the  following  channel -stabi  1  i ty 
related  applications:  (a)  Consolidate  the  channel  model  with 
continuous  sediment  yield  and  bank-stability  models.  Run  the 
consolidated  models  for  a  period  of  a  few  years,  and  predict  the 
size  and  grade  of  channel  needed  to  maintain  a  bank  height-slope 
that  is  stable  for  a  given  stratigraphic  condition.  (b)  Run  the 
consolidated  model  for  a  combination  of  unstable  bank  and  steep 
grade  and  observe  what  combination  of  bed  armoring  and/or  grade- 
control  structures  are  predicted  to  stabilize  the  channel.  (c) 
Select  a  range  of  storm  events  and  use  the  consolidated  model  to 
study  slough  of  bank  material,  and  find  channel  width  and/or 
armoring  coat  that  is  needed  to  prevent  erosion  of  slough  material. 

4.  It  is  recommended  that  data  gathering  efforts  be  continued  to 
provide  an  adequate  base  lor  further  model  development  and 
va 1 idat ion . 
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ADDENDUM  1:  SEDIMENT  CONTINUITY  EQUATION 


Consider  an  unsteady  streamflow  carrying  sediment  down  a  channel 
with  arbitrary  cross-sectional  geometry  and  alignment  (Fig.  1.1).  The 
volume  concentration  of  the  sediment-laden  flow,  c,  may  vary  from  point 
to  point,  and  as  a  result  of  convection,  from  instant  to  instant  at  any 
point.  The  lateral  volume  rate  of  sediment  inflow,  q^,  can  vary  with 
both  space  and  time.  The  continuity  equation  for  an  infinitesimal  unit 
volume  in  the  neighborhood  of  a  fixed  point  is 


3c 

3t 


+ 


=  o. 


(1.1) 


in  which  u^  is  the  instantaneous  local  velocity  component  of  the 
sediment-laden  fluid  along  the  x^-direction .  Repetition  of  the  sub¬ 
script  k  in  a  term  implies  summation  over  the  three  orthogonal  coordi¬ 
nate  directions.  In  a  turbulent  flow,  time  averaging  Eq.  1.1  yields 
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where  the  overbars  indicate  temporal  means,  and  the  primed  terms  repre¬ 
sent  turbulent  fluctuations.  Let  assume  the  mean  cross  sectional  area, 
A,  of  the  channel  divided  in  two  parts.  One  part,  A^ ,  is  occupied  by 
the  sediment  carried  in  suspension,  and  another,  A^,  is  occupied  by  the 
sediment  transported  as  bedload.  It  should  be  noted  herein  that  the 
orientation  of  the  coordinate  system  is  arbitrary  and,  in  general,  A  is 
a  function  of  as  well  as  of  time.  For  convenience,  the  coordinate 
system  will  be  chosen  so  that  the  direction  normal  to  A  coincides  with 
the  streamwise  direction  x  (Fig.  1.1).  The  continuity  equation  for  the 
suspended-load  section  is  obtained  by  integrating  Eq.  1.2  over  A^,  to 
obtain 
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Using  Leibnitz’s  rule: 
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JJ  dAj  -  [c  +  3X  //(CjU  +  cju’J  dAj 
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I(clu  +  C1U’J  9x~t  =  °> 


(1.4) 


in  which  is  the  mean  perimeter  bounding  the  mean  area  A^ . 
dary  condition  for  this  cross  section  is 


The  boun- 


dA 
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9A, 


[ci  dT]  =  lci  (aT  +  11  3xI)]  =  -1'  ci  %  dV 


(l.s) 
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where  q^,  the  volume  rate  of  lateral  sediment  inflow  per  unit  length  of 
,  is  taken  positive  for  inflow.  Expanding  Eq.  1.5  and  time  averaging 
gives 


3Aj  _  9A  __  3A 

["i  aF  +  "i  “  ax  +  w'  9x 


X  (c! +  cjqp  do] .  (1-6) 
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Substituting  this  equation  into  Eq.  1.4  yields 


\l  SI  CjdAj*  SJ  (CjU  +  cju,)dA]  =  /  ( c  j  q^  +  rjqpdOj  (1.7) 


1 


1 
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Analogous  to  the  Reynolds  closure  scheme  for  turbulent  momentum  trans¬ 
fer,  an  eddy  mass  diffusivity  tensor,  e,  is  introduced  such  that 


c:u  =  -e 


3x 


(1.8) 


Now  let 


C]  -  +  c», 


(1.9) 
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and 


u  =  u  +  u* , 


(1.10) 


where  the  tildes  denote  spatial  averages  over  A ^  ,  and  the  asterisks 
denote  deviation  from  the  spatial  averages.  Substitution  of  Eqs .  1.8, 
1.9,  and  1.10  into  Eq.  1.7  gives 
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where 

Qj  =  dA2 

A1 


[  (clq£  +  clq£)  dV 
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(1.11) 


represents  the  volume  rate  of  suspended  sediment  discharge.  The  first 
integral  on  the  right  hand  side  of  Eq.  1.11  denotes  the  rate  of  longitu¬ 
dinal  dispersion  of  sediment.  The  second  integral  can  be  replaced  by 

the  sum  of  the  sediment  inflow  from  runoff  and  tributaries,  q  ,  and  the 

*s 

rate  of  sediment  volume  transfer  across  H-H'  (Fig.  1.1)  from  the  bedload 
zone,  denoted  herein  as  -  Or^/Ot.  Longitudinal  dispersion  in  the  sus¬ 
pended  zone  is  negligible  with  respect  to  vertical  dispersion  and  can  be 
omitted.  Thus,  Eq.  1.11  becomes 


9V~1 

8Qj 

3r 
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Similarly,  integrating  Eq.  1.2  over  A^  yields  the  continuity  equation 
for  the  bedload, 


3  Vi  *2 

3t  +  3x 


J  (?2q2  +  e'q')  do2, 


(1  - 13) 
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where  is  the  volume  rate  of  bedload  discharge,  and  is  the  mean 
perimeter  bounding  .  Two  different  sources  contribute  to  the  integral 
in  the  right  hand  side  of  Eq.  1.13.  One  is  the  rate  of  sediment  volume 
transfer  across  H-H'  from  the  suspended  zone,  and  denoted  Sr^/St.  The 
other  contribution  comes  from  the  rate  of  sediment  exchange  with  the  bed 
surface,  which  can  be  approximated  by 

qb  =  -(l-A)W  ||  ,  (1.14) 


where  A  is  the  bed  porosity,  W  defines  the  active  bed  width  (Fig.  1.1), 
and  z  is  the  local  bed  elevation  above  a  reference  datum.  The  negative 
sign  in  Eq.  1.14  accounts  for  the  fact  that  the  sediment  flux  across 
is  negative  when  sediment  settles  out  of  the  bedload  zone  during  bed 
aggradation  (3z/3t>0).  Replacing  the  integral  in  Eq.  1.13  by  the  sum  of 
S^/Bt  and  q^  gives 


3A  c  3Q  3r  . 

-.-2  2  +  __ 2  =  — 2  -  (i-X)W  — 

3t  3x  3t  U  Al  3t 
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Adding  Eqs .  1.12  and  1.15  results  in 


3Ac  9Qs  n  3z  f3rl 
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where  Ac  =  A^Cj  +  aua  Qs  =  Qj  +  ^2 '  The  tern)  within  brackets 
represents  the  net  rate  of  sediment  volume  transfer  across  H-H' .  This 
term  is  different  from  zero  whenever  a  net  amount  of  sediment  passes 
from  the  bedload  to  the  suspended  load  zone,  or  vice  versa,  as  a  result 
of  bed  scour  or  sediment  deposition.  Letting 


3r 

3t 


and  dropping  the  tildes  and  overbars,  Eq .  1.16  becomes 


3Ac  (  ,  ,  3z  3r  _ 

9t~  +  3x  +  (1*A)  W  3l  +  3t  ~  V 


(1.17) 
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in  which  Qg  is  the  total  volume  sediment  discharge,  and  c  is  the  total 
average  volume  concentration.  Eq.  1.17  is  the  sediment  continuity 
equation  used  in  the  present  model,  and  equally  applies  to  streamflows 
carrying  sediment  mostly  in  suspension  or  as  bedload. 
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ADDENDUM  2:  SUMMARY  OF  TRANSPORT  FORMULAS 
This  Addendum  describes  briefly  the  transport  formulas  used  in  the 
present  model.  Each  formula  is  presented  as  the  originator  intended, 
but  rewritten  in  terms  of  dimensionless  parameters.  Where  the  formula¬ 
tions  required  graphical  solutions  (i.e.,  determination  of  threshold 
conditions  from  Shield's  curve),  analytical  equivalents  (not  shown  here) 
have  been  worked  out  to  facilitate  their  use  in  digital  computation. 
Total  load  formula  of  Yang  (1973): 

*  =  ej*0  Z50  (V/n*)[l0*  *6/S),  (2.1) 

where 


$  =  5.435  -  0.286  log  (w  d^Q/v)-0.457  log  (u^./w)  + 
[1.799  -  0.409  log  (w  d5  /v)  -  0.314  log  (u*/w)] 
log  (VSQ/w  -  VcSQ/w), 


v  J  2.5/[log(u*d5Q/v)  -  0.06]  +  0.66,  0<uJ.d5Q/v<70, 

^  2.05,  u*d50/v  >  70. 

DuBoys-type  bedload  formula  (Graft,  1971): 


(2.2) 

(2.3) 

(2.4) 


where 


<t>  =  a  e2(i  -  e/e  ) 


a-x(S-l)3/Vg3/2d1/2 


(2.5) 

(2.6) 


Empirical  relationships  for  the  sediment  coefficient  x  ate  given  by 
Graft  (1971). 

Bedload  formula  of  Meyer-Peter  and  Muller  (1948): 

<}>  =  8  e3£2  [k3/2  -  ec/em]3/2,  (2.7) 

in  which 

K  =  (V/u*)(fb/8)1/2  (2.8) 

The  friction  factor  associated  with  bed  skin  friction,  f^,  is  obtained 
from  (Vanoni,  1975) 
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(fb)  ’  =  2  log  (4Rb/ks)  +1.14,  — *  >  70, 


(2.9) 

(2.10) 


where  kg  =  dg^,  Np  is  the  flow  Reynolds  number,  and  Rb  is  the  bed  hy¬ 
draulic  radius. 

In  Eqs .  2.1,  2.5,  and  2.7  4>  is  the  dimensionless  volume  transport 
rate,  0^  and  are  the  mobility  number  and  relative  roughness  based  on 
the  d^  grain  size,  and  the  subscript  c  denotes  threshold  conditions. 
These  parameters  are  defined  as 


*  =  CQS/W) / [ (S- 1 ) gd|Q ] ^ , 

(2. 

11) 

ek  =  4/[S-l)gdk], 

(2. 

12) 

and 
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where  S 

is  the  specific  gravity  of  the  bed  material,  and  u, 

j  s  the 

bed 

shear  velocity. 
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ADDENDUM  3:  DESCRIPTION  OF  COMPUTER  PROGRAM 

The  computer  program  used  in  the  simulation  of  the  East  Fork  River 
data  is  listed  below.  The  program  consists  of  a  main  program  and 
several  subroutines.  The  main  program  inputs  the  required  data  to  the 
system,  calls  the  subroutines  according  to  the  computational  scheme,  and 
prints  out  the  calculated  results.  Subroutine  WROUT  is  a  user  supplied 
program  that  routes  water  through  a  channel  reach  for  each  time  step. 
Subroutine  SROUT  performs  sediment  routing  through  a  reach  for  each  time 
step,  and  it  calls  in  turn  subroutines  DUBOYS,  MEYER,  SETVEL,  SHIELD, 
and  YANG.  These  subroutines  compute  potential  carrying  capacities  and 
sediment  transport  parameters. 

In  order  to  minimize  the  need  for  core  memory,  the  program  was 
organized  by  taking  advantage  of  the  absence  of  backwater  conditions  in 
the  present  simulations.  In  the  adopted  organization  the  length  of  the 
entire  channel  reach  is  divided  into  a  few  segments,  or  channel  blocks, 
and  the  time  steps  used  for  the  whole  simulation  period  are  grouped  into 
several  consecutive  time  blocks.  Then,  the  channel  blocks  are  processed 
in  sequence  for  each  time  block.  Input  data  for  each  channel  and  time 
block  are  stored  in  disk  files.  A  list  of  the  important  variables  in 
the  computer  program  is  given  in  the  following  section. 

The  codes  requires  72,741  words  on  a  Mod  Comp  Classic  computer 
system.  This  is  a  16-bit  machine  that  uses  two  words  for  each  single 
precision  variable.  The  execution  time  for  the  East  Fork  River  test  is 
approximately  9  minutes. 

List  of  Fortran  Variables 

Name  Description  Units 

ACCM(IF,  IFA)  Element  of  volume  entrainment  matrix  ft3 

at  the  start  of  the  current  time  step. 

IFA  indicates  the  material  fraction 
exposed  by  the  removal  of  fraction  IF. 

AE  Area  of  flow  cross  section  at  the  ft2 

start  of  the  current  time  step. 
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ARMHT 

BEDELV(IC) 

BEDMAT(IF) 


BEDUP 


BEDWN 


CAP 


CC(IC,  IF) 


CDEP(IC) 


CHEZY 

CHI 

CICIC,  IF) 


Current  thickness  of  active  layer. 

Bed  elevation  of  section  IC  at 
the  end  of  the  current  time  step. 
Vector  used  to  store  the  volume  of 
individual  material  fractions  present 
in  the  bed  layer. 

Change  in  bed  elevation  caused  by 
deposition  during  the  current  time 
step . 

Change  in  bed  elevation  caused  by 
erosion  during  the  current  time  step. 
Volume  concentration  of  individual 
material  fractions  at  transport 
capacity. 

Volume  concentration  of  material 
fraction  IF  passing  through  the 
section  IC  at  the  end  of  the  current 
time  step. 

Coefficient  of  power  formula  relating 
area  of  cross-section  IC  to  water 
depth. 

Chezy's  roughness  coefficient. 
Coefficient  for  the  DuBoys  sediment 
transport  formula. 

Volume  concentration  of  material 
fraction  IF  passing  through  the 
section  IC  at  the  start  of  the 
current  time  step. 


ft 

ft 


ft3 


ft 


ft 


ft3/ft3 


ft3/ft3 


L 

ft  V  sec 


ft5-25/ 


ft3/ft3 
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COEF(J) 


Coefficient  of  the  entrainment 


CONC 

CPER(IC) 

CUP(1T, 

DARM 

DELT 

DELX 

DEPO 

DMM(IF) 

DPTH 

DTS 

EDEP(IC) 


IF) 


frequency  matrix. 

Capacity  concentration  predicteii  by 
transport  formulas. 

Coefficient  of  power  formula  relating 
wetted  perimeter  of  section  IC  to 
flow  cross-sectional  area. 

Volume  concentration  of  inflow  of 
material  fraction  IF  at  the  start 
of  the  current  time  step. 

Size  of  smallest  bed  material  fraction 
that  becomes  immobile  in  the  active 
layer . 

Time  taken  by  a  sediment  characteris¬ 
tic  to  travel  the  distance  DELX. 
Channel  length  increment  used  in  the 
computational  grid. 

Volume  of  fraction  IF  deposited  on 
the  bed  during  the  current  time 
step . 

Representative  size  of  sediment 
fraction  IF. 

Water  depth. 

Size  of  current  time  step. 

Exponent  ol  power  formula  relating 
area  of  cross-section  IC  to  water 
depth 


ppm 


ft3/ft3 


mm 


sec 


ft 


ft3 


mm 


ft 

sec 
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EPER(IC) 


Exponent  of  power  formula  relating 


ERS(IF) 


ERO 


G(IT,  I,  IF) 


GAMA 

GMES(IT,  I,  IF) 


GTOTC(IT) 


GTOTM(IT) 


I 

IF 


INLS 


wetted  perimeter  of  cross-section  IC 
to  flow  cross-sectional  area. 

Volume  of  sediment  fraction  IF 
eroded  from  the  bed  during  the 
current  time  step. 

Parameter  controlling  detachment  of 
bed  material. 

Calculated  volume  discharge  of  sedi¬ 
ment  fraction  IF  passing  through  the 
downstream  end  of  channel  block  I, 
at  the  end  of  the  current  time  step  IT. 
Specific  weight  of  water. 

Measured  inflow  of  sediment  fraction 
IF  to  channel  block  I,  at  the  start 
of  the  current  time  step  IT. 

Calculated  sediment  discharge  passing 
through  the  channel  outlet  at  the 
end  of  the  current  time  step  IT. 
Measured  sediment  discharge  passing 
through  the  channel  outlet  at  the 
start  of  the  current  time  step  IT. 

Index  identifying  channel  block. 

Index  identifying  sediment-size 
fraction . 

Number  of  channel  subreaches  in  a 
channel  block. 


ft3 


ft3/sec 


lbs/ft3 

lbs/sec 


lbs/sec 


lbs/sec 


J.  68 


IT 


Computation  time  step. 


ITCOM 

KIN 

GLAT(IT) 


MEYERP 


N 


NARM 


NFR 

NSEG 

PCBC(IC,  IF) 


PCF ( IC ,  IF) 


PCFF(I ,  IC,  IF) 


Number  of  time  steps  in  simulation 
period . 

Kinematic-wave  routing  parameter. 

Lateral  volume  inflow  of  sediment  ft3/sec 

to  channel  block,  during  the  current 
time  step  IT. 

Coefficient  in  the  Meyer-Peter  and 
Muller  sediment  transport  formula. 

Number  of  time  blocks  the  simula¬ 
tion  period  is  divided  into. 

Integer  identifying  the  smallest 
size  fraction  that  becomes  immobile 
in  the  active  layer. 

Number  of  representative  size 
fractions  used  in  the  simulation. 

Number  of  channel  blocks. 

Percentage  of  material  fraction  IF 
present  in  the  active  bed  layer  of 
cross  section  IC. 

Percent  of  material  finer  than  size 
IF  present  in  the  active  bed  layer 
of  cross  section  IC. 

Initial  percent  of  material  finer 
than  size  IF  present  in  the  active 
bed  layer  of  cross  section  IC  in 
the  channel  block  1 . 
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PCW(IF) 


Percent  of  material  in  transport 
finer  than  size  IF  passing  through 
the  channel  outlet. 


POR(IF) 

Porosity  of  bed  material  fraction  IF. 

ft3/ft3 

Q(IT) 

Computed  water  discharge  at  the  end 

of  the  channel  block  I,  and  at  the 

end  of  the  current  time  step  IT. 

ft3/sec 

QLAT(IT) 

Lateral  water  inflow  to  channel  block, 

during  the  current  time  step  IT. 

ft3/sec 

QMES(IT,  I) 

Measured  water  inflow  to  channel 

block  I,  at  the  start  of  the 

current  time  step  IT. 

ft3/sec 

QUP(IT) 

Upstream  water  inflow  to  channel 

block,  at  the  start  of  the 

current  time  step  IT. 

ft3/sec 

RESCAP 

Residual  transport  capacity  of  an 

individual  material  fraction, 

ft3/ft 

- 

expressed  as  volume  of  dry  sediment 

per  unit  length  of  channel. 

RHB 

Hydraulic  radius. 

ft 

SCAP(IF) 

Potential  transport  capacity  of 

material  fraction  IF,  expressed  as 

volume  of  dry  sediment  per  unit 

length  of  channel. 

ft3/ft 

SLN 

Length  of  channel  block. 

ft 

SLOPE (IC) 

Channel  bed  slope  at  section  IC. 

ft/ft 
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SHU 

Average  value  of  the  kinematic 

viscosity  of  water  for  the  simu¬ 
lation  period. 

f t2/sec 

SPGR 

Specific  gravity  of  sediment  material. 

- 

SUM 

Summation  term  in  Eq.  11. 

- 

TAO 

Average  unit  tractive  force. 

lbs/ft2 

TBM 

Total  volume  of  bed  materia] 

contained  in  the  active  layer  per 

unit  length  of  channel. 

ft3 

TC 

Critical  unit  tractive  force. 

lbs/ft2 

TCA(IF) 

Sum  of  all  the  elements  in  row  IF 

of  the  volume  entrainment  matrix. 

ft3 

TEMP 

Average  water  temperature  for  the 

simulation  period. 

Fahrenhei 

TGIN(IF) 

Vector  used  to  store  the  measured 

volumes  of  all  fractions  that  enter 

the  channel  during  the  simulation 

period. 

ft3 

TGM(IF) 

Vector  used  to  store  the  measured 

volumes  of  all  fractions  that  leave 

the  channel  during  the  simulation 

period . 

ft3 

TGMES 

Total  sediment  yield  measured  at  the 

channe 1  outlet. 

lbs 

TGO(IF) 

Vector  used  to  store  the  computed 

yield  of  individual  fractions. 

lbs 
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TGOUT 


Total  sediment  yield  computed  at 
the  channel  outlet. 


lbs 


UST 

VEL 

VS(IF) 

WATMAT(IF) 

WEP 

WIDTH (IC) 

XSI(IC) 


Bed  shear  velocity. 

Average  velocity  of  flow. 

Settling  velocity  in  quiescent  water 
for  material  fraction  IF. 

Vector  used  to  store  the  volumes  of 
all  the  material  fractions  being 
carried  by  the  flow  during  the 
current  time  step. 

Wetted  perimeter  of  channel  cross 
section. 

Vector  used  to  store  the  active  bed 
width  of  all  channel  cross  sections. 
Distance  of  cross-section  IC  to 
upstream  boundary  of  channel  block 
containing  IC. 


ft/sec 
ft/ sec 
ft/sec 

ft3 


ft 

ft 

ft 


1 


o  n 


C  LIST  OF  COMPUTER  PROGRAM  USED  IN  THE  EAST  FORK  RIVER  TEST 
C  fr************************************************* ******* 

C 

DIMENSION  T ITLE (20) , QMES (264 ,3) ,QOUT(30  0  > ,GMES(264 ,10,3), 

4PCFF  <3,6,10) , PCF  <21 ,10) ,PCFR (21 ,10), GTOTM  <  264 ) , CTQTC ( 264 ) , 

&QS (21 ) ,INT(21) ,GS(21) ,CII <3,6,10) , BEDLVL (3,6) ,TGO( 10  > ,TGM( 10)  , 
6TGIN ( 1 0 ) 

COMMON  /ROUT/  I , IT , SLN , CHEZY , DTS , I TCOM , INL , QBASE , 

&QUP (300),QI(S0,2) , XI (SO ) ,KSI (SO )  ,0(300 ,3) 

COMMON  /WROUT/  SLOP , QLAT ( 30 0 ) , QC ( SO , 2 ) , XC ( SO > , KSC < SO ) 

COMMON  /SROUT/  SP GR , CAMA , SNU , NFR , I NLS , WIDTH ( 1 0 ), SLOPE ( 1 0  )  , 

6CPER (10) , EPER(iO) ,CDEP (10 ) ,EDEP( 10) ,GLAT (300 , 10) ,XSI( 10) , 

&XSC< 10 , 10) , BEDELV ( 10),CI(10,10>,CC(10,10),PCBI(10,10) , PCBC( 10 , 10  >  , 
6CUP (300,10) , DMM (10)  ,COEF( 10) , US ( 10) , POR  < 1 0 ) ,G(264 ,3 ,10  > ,PCW< 10 ) , 
6ER0 , CHI , CEL 
C 

C  DATA  INPUT 
C  GENERAL  INFORMATION 

READ (4,406)  NSEG , DTS , I TCOM ,NFR 
C  PHYSICAL  PROPERTIES 

READ (4,407)  TEMP , GAMA , SPGR  , SNU 
READ (4,408)  ( DMM ( IF > , IF  =  1 , NFR  ) 

C  WATER  AND  SEDIMENT  ROUTING  PARAMETERS 
READ  <4,40V)ERO,CHI,CEL, CHEZY 
C  INITIAL  BED  MATERIAL  SIZE  DISTRIBUTION 

READ( 4 , 40S )  ( (PCFF( 1 , IC, IF >  ,  I F  -  1 , 1  0  >  ,  I  C=  1 , 3  > 

READ (4 ,405)  ( ( P CFF < 2 , I C , IF ) , IF=1 , 1 0 ) , I C= 1 , 4 ) 

READ (4,405)  ( (PCFF (3, IC, IF), IF- 1,10) , IC-1 ,6) 

C  COEFFICIENTS  OF  ENTRAINMENT  FREQUENCY  MATRIX 
COEF ( 1 ) =1 .0 
COEF ( 2 ) =1 . 0 
DO  280  J- 3 , NFR 
280  C0EF(J>-2. OfcCOEF  <  J- 1  ) 

SUBW=  ( SPGR-1  .  0  )  i'GAMA 
C  POROSITY  AND  SETTLING  VELOCITY 
DO  88  1=1, NFR 

POR  < I )  =  1 .-0 .245-0 . 0864/(0 . S*DMM(I >>**0.21 
D=DMM( I) 

CALL  SETVEL ( D ,W ) 

VS<  I )=U/30 . 48 
88  CONTINUE 

BEGINNING  OF  TIME-BLOCK  LOOP 
DO  999  N=i , 4 

C  MEASURED  WATER  DISCHARGE  AT  SECTIONS  B-l ,  B-5M,  AND  B-17 
DO  351  1=1,3 

READ (3, 333)  <  QMES (IT, 1 ) , I T  =  1 , 1 TCOM) 

351  CONTINUE 

C  MEASURED  SEDIMENT  DISCHARGE  AT  SECTIONS  B-l,  B-5M,  AND  B-17 
DO  350  11=1,3 
DO  888  IT= 1  ,  ITCOM 
QLOG=ALOGi  0  ( QMES  ( IT ,  1. 1  )  ) 

IF(Il.EQ.i)  EXPONT=3.05785*QLOG-9  33462 
IF (II  . EQ.2)  EXPONT=l . 92638*QL0G-3 . 294 09 
IF( II . EQ. 3)  EXPONT  =  1 . <16664*G)L0G-S  0  3701 
GTOTM( IT )=1 0 . 0**EXPONT 
888  CONTINUE 

READ( 4 ,401)  NINT 
WR ITE ( 5 , 40 1 )  NINT 

READ (4 ,402)  <  INT  ( I )  ,  G|S<  I  >  ,  GS  ( I  >  ,  (  PCF  ( I  ,  J  )  ,  J  =  1  ,NFR  )  ,1  =  1  ,NINT> 
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WRITE <5, 402)  <INT<I) ,QS(I ) ,GS(I ) , < PCF ( I , J ) , J  =  1 ,NFR ) ,1  =  1 ,NINT) 
12=11 

DO  309  1=1 ,NINT 
PCP  =  0 . 0 

DO  301  J  =  1  ,NFR 
PC=PCF  ( I , J ) -PCP 
PCFR ( I , J )=PC/1 00 . 0 
PCP=PCF ( I  ,  J) 

301  CONTINUE 
309  CONTINUE 

DO  302  IF=i,NFR 
DO  303  1=1 ,NINT 
IF(I.GT.l)  GO  TO  304 
INT1  =  INT  ( 1 ) 

DO  305  J=i , INT i 

305  GMES( J, IF , II )=PCFR  (l,IF)*GTOTM(J) 

KC=INT  < 1 ) 

GO  TO  303 

304  IF(I  EQ.NINT)  GO  TO  306 
IB=KC+i 
IE  =  INT  < I ) 

GD= (PCFR ( I , IF ) -PCFR ( 1-1 , IF ) ) /FLOAT ( IE-KC) 

DO  307  J=IB , IE 

307  GMES<  J  , IF , I 1 )  =  (PCFR  ( I -1 , IF  )+GD#FLOAT  (  J-KC) )  #GTOTM(  J  ) 

KC=INT  < I ) 

GO  TO  303 

306  CONTINUE 
INT2=INT (NINT ) 

DO  308  J=INT2, ITCOM 

308  GMES( J ,  IF  ,  1 1 )=PCFR (NINT , IF)*GTOTM( J ) 

303  CONTINUE 

302  CONTINUE 
350  CONTINUE 
C 

C  BEGINNING  OF  CHANNEL-BLOCK  LOOP 
DO  201  1=1 , NSEG 
INL=4  0 

READ (1,403)  INLS , SLN 
WRITE (5, 403)  INLS, SLN 
C  GEOMETRIC  INPUT  FOR  CHANNEL  BLOCK 

READ(1 ,404)  ( XSI ( IC ) , SLOPE < IC ), WIDTH ( I C >, CPER ( IC ), EPER ( IC ) , 
&CDEP ( IC ) ,EDEP ( IC) ,IC=i ,INLS) 

WRITE (5, 40 4)  (XSI(IC) ,SLOPE(IC) , WIDTH ( IC ), CPER ( IC ), EPER ( IC > , 
&CDEP ( IC ) , EDEP ( I C ) , I C= 1 , INLS  > 

DO  801  IC=1 , INLS 
PCP=0 . 0 

DO  802  IF=1 , NFR 
PCF(IC,IF)=PCFF(I, IC,  IF) 

PCBI(IC,IF) =PCF (IC,IF)-PCP 
PCP=PCF ( IC , IF ) 

802  CONTINUE 
801  CONTINUE 
C 

C  INITIAL  AND  BOUNDARY  CONDITIONS  FOR  CHANNEL  BLOCK 
SLOP=0 . 0 

DO  11  IT=i, ITCOM 
IF(I-2)3,4,S 

3  QUP (IT) =QMES ( I T , i ) 

GO  TO  6 

4  QUP ( IT) =QMES( IT , 2) +Q ( IT ,1) 
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GO  TO  6 

5  QUP  <IT)=Q(IT,2> 

6  QLAT ( IT)  =  0  .  0 
DO  11  IF=i,NFR 
IF<I-2>7,8,9 

7  CUP  <IT,IF)~GMES<IT,IF,1)/<QUP(IT)*GAMA*SPGR> 

GO  TO  10 

8  CUP (IT,IF)=(G(IT,1,IF> +GMES1 IT , IF , 2  > / <GAMA#SPGR ) > /QUP  < IT ) 

GO  TO  10 

9  CUP  <IT,IF)=G<IT,2,IF)/QUP<IT> 

10  GL.AT  (IT, IF) =0.0 

11  CONTINUE 

DO  13  IC  =  i , INLS 
BEDELV  ( IC )  =BEDLVL  ( I  ,  I C ) 

DO  12  IF=i , NFR 

Cl ( IC , IF  > =CI I ( I > IC , I F  > 

IF(CI(IC,IF>  .EQ. 0 . 0)CI(IC,IF)=CUP(1,IF> 

12  PCBC(IC,IF)=PCBI(IC,IF) 

SLOP=SLOP +SLOPE < IC) 

13  CONTINUE 
SLOP=SLOP/FLOAT< INLS) 

QBASE=QUP ( 1  ) 

DO  300  IC=i  ,  INL 
QI(IC,1)=QUP<1) 

QI  <  IC >  2)  =Ql)P  <  1 ) 

XI  <  IC  >  =Sl.N#FLOAT  <  IC -1  > /FLOAT  (INL) 

300  KSI<IC)---0 

C  LIST  INITIAL  BED  ELEVATION  AND  SIZE  COMPOSITION 
URITE(S,50i)  I 

UR ITE< 5,502)  < IC , < PCF < IC , IF ) , IF= 1 , NFR >, BEDELV ( IC )> I C=1 , INLS ) 
C 

C  ROUTE  UATER  AND  SEDIMENT  THROUGH  CHANNEL  BLOCK  DURING 
C  CURRENT  TIME  STEP 

DO  101  IT  =  1 , ITCOM 
C 

C  UATER  ROUTING 
CALL  UR0UT2 
C 

C  SEDIMENT  ROUTING 
CALL  SR0UT2 
C 

IF < INL . EQ . 0 )  GO  TO  129 

C  RESET  INITIAL  CONDITIONS  FOR  FLOW  CALCULATIONS 
DO  132  J  =  1  , 1  Nl. 

XI(J>=XC( J) 

QI ( J , 1 ) =QC ( J , 1 ) 

QI  <  J ,  2 )  =  QC  (  J ,  2 ) 

KSI <  J ) =KSC  <  J) 

132  CONTINUE 
12.9  CONTINUE 

C  RESET  INITIAL  CONDITIONS  FOR  SEDIMENT  CALCULATIONS 
DO  270  IF=1 , NFR 
DO  271  IC=i  ,  INLS 
IFUC.EQ.  1)  GO  TO  272 
CI<IC,1F)=CC(IC,IF) 

GO  TO  271 

272  CI(IC,JF)  =-CUP  <  IT  ,  IF  ) 

271  CONTINUE 
270  CONTINUE 

IF1IT.NE.170)  GO  TO  101 
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WR ITE ( 5 ,S0  3 )  IT 
DO  900  IC=i  , INLS 
PCF(IC,i)=PCBC(IC,l> 

DO  901  IF=2,NFR 

901  PCF(IC,IF)=PCF(IC,IF-1) +PCBC ( IC , IF ) 

900  CONTINUE 

C  LIST  BED  ELEVATION  AND  SIZE  COMPOSITION  AT  TIME  STEP  170 

WRITE (5,S02>  (IC, (PCFtIC, IF) , IF=i , NFR >, BEDELV < IC) , IC=i ,INLS) 
101  CONTINUE 
C  DISCRETIZED  OUTFLOW 
QP=QBASE 

DO  315  IT=1 , ITCOM 
QCC=Q(IT,1) 

Q(IT,I)=(QP+QCC)/2. 0 
QP=QCC 

31S  CONTINUE 
C 

WRITE  <  5  f 504 ) 

DO  911  IC=1 , INLS 
PCF(IC,i)=PCBC(IC,l> 

DO  912  IF  =  2 , NFR 

912  PCF(IC,IF)=PCF(IC,IF-1 >+PCBC( IC , IF ) 

911  CONTINUE 

C  LIST  FINAL  BED  ELEVATION  AND  SIZE  COMPOSITION 

WRITE( 5 , 5021  (IC,  (PCFdC,  IF )  ,  IF=1 ,  NFR  > ,  BEDELV  (IC ) ,  IC=1 ,  INLS  > 
C 

DO  913  IC=i  ,  INLS 
BEDLVL  (1,10=  BEDELV  ( I C ) 

DO  914  IF=1,NFR 
PCFF(I,IC,IF)=PCF(IC, IF) 

CII  (I  ,IC,IF)=CI(IC,IF) 

914  CONTINUE 

913  CONTINUE 
201  CONTINUE 

C  END  OF  CHANNEL-BLOCK  LOOP 
REWIND  1 

C  COMPUTE  HYDROGRAPH,  SEDIMENTGRAPH ,  AND  SEDIMENT  YIELD  AT 
C  SECTION  B-17 
TGOUT=0  0 
TGMES=0  0 
DO  299  IF  =  1  ,NFR 
TGO ( I F ) =0  0 
TGM ( IF ) =0 . 0 
TGIN< IF ) =0 . 0 
299  CONTINUE 

WRITE (5, 505) 

DO  41  IT=1, ITCOM 
QOUT  <IT)=Q(IT , NSEG ) 

GTOTC ( IT ) =0 . 0 
DO  42  IF=i , NFR 

G< IT , NSEG , IF)=G( IT ,NSEG , IF) #GAMA#SPGR 
GTOTC(IT)=GTOTC(IT)+G(IT,NSEG,IF> 
TGO(IF)=TGO(IF)+G(IT,NSEG,IF)*DTS 
TGM<IF)=TGM(IF) +GMES (IT, IF, 3) *DTS 

TGIN( IF)=TGIN(IF)+(GMES(IT ,IF , 1 ) +GMES ( IT , IF , 2 > ) *DTS 
42  CONTINUE 

TGOUT=TGOUT+GTOTC(IT)*DTS 
TGMES=TGMES+GTOTM ( IT ) *DTS 

C  COMPUTE  AND  LIST  CHANGES  IN  SIZE  COMPOSITION  OF  SEDIMENT 
C  LOAD  AT  B-l7  DURING  CURRENT  TIME  BLOCK 


990 


99  i 


992 


43 


902 

41 

C 


999 

C 

C 

333 

401 

402 

403 

404 

405 

406 

407 

408 

409 

501 

502 

503 

504 

505 

506 

509 

510 


511 


IF(N.NE.i)  GO  TO  990 

IF(IT .EQ.213 .OR . IT.EQ .214 . OR . IT .EQ.264)  GO  TO  43 
GO  TO  41 

IF(N  . NE . 2  >  GO  TO  991 

IF<  IT .EQ. 3  OR  IT . EQ . 4  OR . IT  EQ . 54  OR . IT  EQ .99)  GO  TO  43 
IFdT.EQ.100  OR  IT.EQ.144.0R.IT.EQ.148.0R.IT  EQ.149)  GO  TO  43 
IF(IT. EQ. 1 89. OR  IT.EQ.191.0R.IT  EQ. 193. OR. IT. EQ. 239)  GO  TO  43 
IFCIT . EQ  240  OR . IT . EQ . 241 . OR . IT . EQ . 244 . OR  .  IT . EQ .247)  GO  TO  43 
GO  TO  41 

IF < N . NE . 3 )  GO  TO  992 

IF(IT  EQ  21  OR . IT . EQ  22  OR . IT .EQ.24  OR . IT  EQ.69)  GO  TO  43 
IF < IT . EQ . 70 . OR  IT . EQ . 72 . OR . IT . EQ . 1 19 . OR . IT . EQ . 120 >  GO  TO  43 
l!'<IT. EQ. 165. OR  IT  EQ. 167. OR. IT.EQ. 212. OR. IT. EQ.213)  GO  TO  43 
IF(IT.EQ. 262. OR. IT.EQ. 263. OR  IT. EQ.264)  GO  TO  43 


GO  TO  41 

IFdT.EQ  57. OR  .  IT  .  EQ  100  OR  IT  .EQ  144  OR  .  IT  EQ  185)  GO  TO  43 
IF(IT .EQ . 186  OR  IT . EQ  192  OR . IT  EQ . 237  OR . IT . EQ . 238)  GO  TO  43 
GO  TO  41 
CONTINUE 

PCW(1)=G(IT,NSEG,1) /GTOTC (IT) *100  .  0 
DO  902  IF  =  2 , NFR 

PCM  < IF ) =PCW  <IF-i)+G<IT, NSEG , I F ) /GTOTC  < IT ) #1 0  0 . 0 
WRITE <5, 506)  N , I T , ( PCW ( IF ) , IF  =  1 , NFR ) 

CONTINUE 

LIST  SEDIMENT  YIELD,  HYDROGRAPH,  AND  SEDIMENTGR APH  AT  8-17 
WR ITE  <  5 , 509 )  TGMES , TGOUT 
WRITE <5, 510) 

UR ITE <5, 511)  < IF,TGO< IF) , TGM(IF) ,TGIN(IF) , IF=i  , NFR  ) 

WRITE (5,512) 

WRITE<S,513)  <1, (QMESII, J) , J=i ,3),  QOUT ( I ) , GTDTM ( I ) , GTOTC ( I ) , 
41=1 , ITCOM ) 

CONTINUE 

END  OF  TIME-BLOCK  LOOP 


FORMAT ( 1 0F7 .2) 

FORMAT  < 14 ) 

FORMAT <14 ,F6 . 1 ,F8  3, 10F6 . 2) 

FORMAT  < 14 , Fi 0 .2) 

FORMAT  <  F7 . 1 , F3 . 5 , F7 . 3 , F7 . 3 , 3F6 . 3 ) 

FORMAT <10F7. 2) 

FORMAT (14, F7 .1,14,14) 

FORMAT(3F7.2,FiO .7) 

FORMAT  < 1 0F7 . 2 ) 

FORMAT  <  4F 1 0 .5) 

F0RMAT(//15X, 'INITIAL  SIZE  DISTRIBUTION  OF  BED  MATERIAL  AND', 

6'  BED  ELEVATION  FOR  EACH  SECTION  OF  CHANNEL  BLOCK ',14/) 

FORMAT  < 15X , 1 4 , 1 0F8 . 3 , F 15 . 7  ,  '  FEET') 

F0RMAK/15X, 'SIZE  DISTRIBUTION  OF  BED  MATERIAL  AND  BED  ELEVATION 
4'  FOR  EACH  SECTION  OF  CHANNEL  BLOCK  AT  TIME  STEP', 15/) 
F0RMAT(/15X, 'FINAL  SIZE  DISTRIBUTION  OF  BED  MATERIAL  AND  BED', 

4'  ELEVATION  FOR  EACH  SECTION  OF  CHANNEL  BLOCK'/) 

F0RMAT<//15X, 'VARIATION  OF  LOAD  SIZE  DISTRIBUTION  AT  B-17  ', 

4 '  DURING  CURRENT  TIME.  BLOCK'/) 

F0RMAT(5X , 'TIME  BLOCK ', 12 , 3X> ' TIME  STEP ' , 1 4 , SX , 1 0F8 . 3 ) 
F0RMAT<//15X , 'TOTAL  MEASURED  SEDIMENT  YIELD  =',F20.S,'  LBS.'/ISX 
4 ' TOTAL  COMPUTED  SEDIMENT  YIELD  =',F20  5,'  LBS.'//) 

F0RMAT(//13X, 'SEDIMENT  COMPUTED  MEASURED 

4,'  MEASURED' /13X, 'TRACTION  YIELD(LBS)  YIELD', 

4 ' (LBS)  I NFL  OW (LBS ) ' / ) 

FORMAT ( 15X, 13, F20 . S , F20 . 5 , F2 0  5 > 
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512  F0RMAT<//S4X, 'MEAS. '/16X, 'TIME  COMPUTED  FLOW(CFS)  AT 

& ,  ' FLOW  SED.  LOAD(L.BS)  AT  B-17'/16X,  'STEP  B-l  B-SM', 

t.'  B-17  B-17  MEAS.  COMP.'//) 

513  FORMAT  < 15X , 14 , 6F1 0 . 3 ) 

520  FORMAT(2X,F20 .5/2X,E15.7) 

END  FILE  5 
STOP 
END 
C 

SUBROUTINE  WR0UT2 
C 

C  THIS  SUBROUTINE  ROUTES  WATER  THROUGH  A  CHANNEL  BLOCK  USING  THE 
C  KINEMATIC  WAVE  SCHEME  DEVELOPED  BY  BORAH  ET  AL<1980).  A  LIST  OF  THE 
C  IMPORTANT  VARIABLES  USED  IN  THIS  SUBROUTINE  IS  GIVEN  BELOW: 

C  - 

C  NAME  DESCRIPTION  UNITS 

C  - 

C  IC 

C 

C  INL 
C 

C  IS 

C 
C 

C  KSC(IC) 

C 
C 
C 

C  KSI(IC) 

C 
C 
C 

C  Q ( IT , I ) 

C 
C 

C  QC(IC,L) 

C 
C 

C  Q I  ( I C  ,  L. ) 

C 
C 

C  XC(IC) 

C 
C 

C  XI(IC) 

c 
c 

c  - 

c 

COMMON  /ROUT/  I  ,  IT  , SLN  , CHEZY  , DTS  ,  ITC.OM  ,  I NL  , QBASE  , 

&QUP(300) ,QI(50,2) , XI < SO  ) , KSI ( 50 ) , Q < 30 0 , 3 ) 

COMMON  /WROUT/  SLOP , QLAT < 300 > , QC < 50 , 2 > , XC ( 50 > , KSC < 50 > 

EXP-1  S 
BET=1 . 0/EXP 
EXP  1 =EXP  + 1 . 0 
EXMl=EXP-i . 0 
KIN=CHEZY*SLOP**0 .5 
TF.RM=EXP*KIN*DTS 
QL-QLAT  < IT ) 

QU=QUP (IT) 


PARAMETER  USED  TO  LABEL  INDIVIDUAL  WAVES. 

NUMBER  OF  NODE  POINTS  IN  THE  CHANNEL  BLOCK 

COUNTS  NUMBER  OF  SHOCK  WAVES  OCCURRING  AT  THE  END 
OF  THE  CURRENT  TIME  STEP  . 

FLAG  USED  TO  CHARACTERIZE  THE  WAVE  IC  OCCURRING  IN 
THE  CHANNEL  BLOCK,  AT  THE  END  OF  THE  CURRENT  TIME 
STEP.  KSC=0  FOR  CHARACTERISTICS,  KSC>0  FOR  SHOCKS. 

FLAG  USED  TO  CHARACTERIZE  THE  WAVE  IC  OCCURRING  IN 
THE  CHANNEL  BLOCK,  AT  THE  START  OF  THE  CURRENT  TIME 
STEP.  KSI =0  FOR  CHARACTERISTICS,  KSI>0  FOR  SHOCKS. 

COMPUTED  WATER  OUTFLOW  FROM  THE  CHANNEL  BLOCK  I  ,  AT  C.FS 
THE  END  OF  THE  CURRENT  TIME  STEP  IT. 

FLOW  DISCHARGE  AHEAD  <L=i>  OR  BEHIND  <L=2)  OF  THE  CFS 
SHOCK  IC,  AT  THE  END  OF  THE  CURRENT  TIME  STEP. 

FLOW  DISCHARGE  AHEAD  (L=l)  OR  BEHIND  <L=2>  OF  THE  CFS 
SHOCK  IC,  AT  THE  START  OF  THE  CURRENT  TIME  STEP. 

DISTANCE  OF  WAVE  FRONT  IC  TO  UPSTREAM  END  OF  THE  FT 

CHANNEL  BLOCK,  AT  THE  END  OF  THE  CURRENT  TIME  STEP. 

DISTANCE.  OF  WAVE  FRONT  IC  TO  UPSTREAM  END  OF  THE  FT 

wHANNEL  BLOCK,  AT  THE  START  OF  THE  CURRENT  TIME 
STEP  . 


IMQL.EQ.  0  .  0  ANP.RU.FR.  0  .  0)  GO  TO  130 
C  PROJECT  ALL  CHARACTERISTICS  TO  NEW  TIME  LEVEL 
AC=(QU/KIN)**BET+QL*DTS/2  0 
QC(i,l)=KIN*AC**EXP 
IF  <  QL . EQ  0 . 0 )  GO  TO  102 
XC<l)=<QC<i,i)-QU>/QL 
GO  TO  103 

102  XC<i>=TERM*AC**EXMi/2  0 

103  ICST  =  1 
KSC ( 1 >=0 
GO  TO  131 

130  ICST=0 

131  IF < INL . EQ . 0 )  GO  TO  ISO 
IS=0 

DO  104  IC=1 , INL 
IA=IC+ICST-IS 

IF(KSKIC)  -EQ.0  OR  QI < IC , 1 )  . GE.QI(IC,2>>  GO  TO  10S 
C  PROPAGATION  OF  SHOCK  WAVE 
AA=(QI(lC,i)/KIN)**BET 
AB- <  QI <IC,2)/KIN)**BET 
AAF=AA+QL*PTS 
ABF  =  AB+QL.*DTS 
QC<  IA,1  )=KIN#AAF**EXP 
QC< IA,2)=KIN*ABF**EXP 
I F ( QL . EQ . 0 .0)  GO  TO  10B 
PROD=ALP/<EXPl*< AB-AA)#QL> 

XC< IA)=XI < IC>+PROD*<ABF**EXPl-AAF*#EXPl-AB**EXPi+AA**EXPi) 
GO  TO  107 

108  XC<IA)=XI<IC)  +  <QI(IC,2)-G)I(IC,1 >  > *DTS/ ( AB-AA > 

GO  TO  107 

C  PROPAGATION  OF  CHARACTERISTIC  WAVE 

105  KSI<IC)=0 

AC=(QI  (IC  ,  1>/KTN>**BET+QL*DTS 
QC( IA, 1 )=KIN*AC*#EXP 
IF  <  QL . EQ . 0 . 0 )  GO  TO  106 
XC<IA)=XI(IC)  +  (QC(IA, 1)  -  QI(IC,1 ))/QL 
GO  TO  107 

106  XC< IA)=XI < IC)+TERM*AC**EXM1 
C  CHECK  FOR  NEW  SHOCK  FORMATION 

107  IF<KSI<IC) GT.0)  GO  TO  10? 

IF < I A . EQ . 1 >  GO  TO  140 
IF<XC(IA) .LE.XC(IA-l) )  GO  TO  110 

C  NO  SHOCK  IS  FORMED 
140  KSC ( I A ) =0 
GO  TO  104 

C  SHOCK  IS  FORMED 
110  IA= I A-l 

IF  <  K  SC ( I A ) . GT . 0 )  CO  TO  112 
C  SHOCK  IS  FORMED  BY  TWO  CHARACTERISTIC  WAVES 
XC<  IA)  =  (XC.(1A)+XC(  IA+i  )  )/2 . 0 
QC<IA, 1 >=QC< IA+i , 1 ) 

QC(IA,2)=QC(IA,1) 

IS--=IS+i 
K  SC  <  I A  )  =  1 
GO  TO  104 

C  THE  CHARACTERISTIC  WAVE  JOINS  THE  SHOCK  AHEAD  OE  THE  FRONT 
112  XC  <  I A  )  =XC  (  I A ) 

QC(IA,1 )=QC< IA  + 1,1) 

IS=IS+i 
K  SC< I A ) =1 1 
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GO  TO  104 

C  CHECK  IF  THE  PROPAGATING  SHOCK  IS  INTERSECTED  BY  ANY  OTHER  SHOCK 
C  OR  CHARACTERISTIC  WAVE 
109  IF  < I A  EQ . 1 )  GO  TO  141 

IF(XC<IA>  LE.XC(IA-D)  CO  TO  113 
C  THE  PROPAGATING  SHOCK  IS  NOT  INTERSECTED 
141  KSC  <  I A  )  =  1 

GO  TO  104 

C  THE  PROPAGATING  SHOCK  IS  INTERSECTED 

113  IA=I A-i 

IF  <KSC< IA> . GT . 0 )  GO  TO  114 

C  THE  PROPAGATING  SHOCK  IS  JOINED  BY  A  CHARACTERISTIC  WAVE 
C  BEHIND  THE  FRONT 

XC< IA)=XC( IA+1 ) 

OC< IA ,2)=QC< I A  ,  1  ) 

QC< IA, 1 )~QC< IA+1 , 1 ) 

ISdS+1 
KSC ( I A  >  ~ 1 1 
GO  TO  104 

C  NEW  SHOCK  IS  FORMED  BY  TWO  INTERSECTING  SHOCKS 

114  XC(IA)=<XC< IAI+XCC IA+1 ) > /2  0 
QC<IA,i)=QC<IA+i  ,1) 

QC(IA,2)=QC(IA,2) 

IS=IS+1 
KSC  < I A ) =2 
104  CONTINUE 

C  COMPUTE  OUTFLOW  FROM  CHANNEL  BLOCK  DURING  CURRENT  TIME  STEP 
ISO  CONTINUE 
IAD=0 

IF(INL.EQ.O)  I A=ICST 
IF<  IA  .  Et) .  0  )  GO  TO  122 
XB=0 . 0 
QB=QIJ 

DO  115  J  =  1  ,  IA 

IF(XC(IA)  .GE.SI.N)  GO  TO  117 
XB---XC(IA) 

QB=QC ( I A , 1 ) 

IF<KSC(IA) .  GT  .  0  )  QB=(QC(IA,l)+QC<IA,2>>/2. 0 

117  IC=IA-J+1 

IF(XC<IC)  LT.SLN)  GO  TO  115 

IAD=IAD+1 

XA=-XC(IC) 

QA-QC (10,1 ) 

IF  <  KSC  < IC )  . GT . 0 )  QA=  <  QC ( IC, 1)+QC<IC,2> )/2. 0 

115  CONTINUE 

IF ( I AD . GT  0 )  GO  TO  118 
AIMQBASE/KIN)**BET 

IF< IT .GT . 1 )  AI=(Q(IT-1 , I ) /K IN ) **BET 

AC=AI+QL*DTS 

QA=KIN*AC**EXP 

XA=TERM*AC*#EXM1 +SLN 

IF  <QL  GT  .  0  .  0)  X  A=  <  QA-Q  ( IT  - 1  ,  I  )  )/Qt„+SLN 

118  Q(IT,I)~QB+<  QA-QB ) *  <  SLN-XB  > / ( XA-XB  > 

INL=IA-I AD 

GO  TO  123 

122  Q( IT  > I  )  =  0  .  0 
INL  =  0 

123  RETURN 
END 

C 
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SUBROUTINE  SR0UT2 


THIS  SUBROUTINE  PERFORMS  SEDIMENT  CALCULATIONS  AT  ALL  NODE  POINTS 
IN  THE  CHANNEL  BLOCK  DURING  CURRENT  TIME  STEP 

DIMENSION  WATMAT ( 1 0  > , BEDM  AT ( 1 0  )  ,SCAP (10>,ERS(10)  , ACCh< i 0 , 10  )  , 

4>TCA(  10  > 

COMMON  /ROUT/  I , IT , SIN , CHEZY , DTS , I TCOM , I NL , QBASE , 

&QUP (300) ,QI(50,2),XI(SO),KSI(SO),Q(300,3) 

COMMON  /SR OUT/  SPGR , GAMA , SNU , NER , I NLS , WIDTH < i 0 ) , SLOPE ( 1 0  )  , 
iCPER < 10 ) ,EPER  < 10  > ,CDEP  < 10 ) , EDEP ( 1 0  > , GLAT  <  3  0 0 , 1 0  ) , XSI ( 1 0  )  , 
iXSCCl 0 , 10 ) , BEDELV (10) ,CI(i0(10),CC (10,10), PCB I<10,10) ,PCBC( 10 , 10) , 
&CUP(300 ,10) , DMM( 10 ) ,COEF ( 10 ) , OS ( 1 0 ) , POR ( 1 0 ) ,G(264,3,10) ,PCU( 10  > , 
4ER0  ,CHI , CEL 

DO  204  IC=i  ,  INI.S 

INTERPOLATE  FLOW  VALUES  AT  THE  NODE  POINTS 
IFdC.EQ.l)  GO  TO  20S 
IF(XSI(IC>  L.T  X I ( 1 ) )  GO  TO  206 
DO  207  J  =  1  ,100 

IF(XSI ( IC) . GE ,XI( J)  AND . XSI ( IC) . LT . XI ( J+i) )  GO  TO  208 
GO  TO  207 

208  Ql=QT ( J  ,  1 ) 

IF (KSI(J)  GT.0)  Ql=(QI(J,i)+QI(J,2))/2.0 
Q2  =  QI  <H  1,1) 

IF  <  KSI ( J  + 1 )  .  GT  0  )  Q2=(QI(J  +  1 ,1 )+QI(J+2,2) )/2. 0 
QP  =  Ql  +  (Q2-Qi)*<XSIdC)-XI(  J) ) / ( XI ( J  +  l ) -XI ( J )  ) 

GO  TO  209 
207  CONTINUE 
GO  TO  209 

206  IFdT.EQ.  1)  GO  TO  210 
Q1=QUP (IT-1  ) 

GO  TO  211 

210  Qi =QI ( IC , 1 ) 

211  Q2=QI(l,i) 

IF (KSI ( 1 ) . GT . 0  >  02- ( QI ( 1 , 1 ) +QI ( 1 , 2 ) > /2 . 0 
RP -QI +  <  Q2-Q 1 ) *XSI < IC ) /XI ( 1 ) 

GO  TO  209 

20S  IFdT.EQ. 1)  GO  TO  212 
G!P=G)UP  (IT-1  ) 

GO  TO  209 

212  QP  =  RI  ( IC  ,  1  ) 

209  CONTINUE 
C 

C  UPDATE  THE  HYDRAULIC  PARAMETERS  AT  THE  NODE  POINTS 
SLP=SL  OPE ( 1C) 

UH)TH=UIDTH(  IC) 

CPR-CPER  (JO 
F-PR-EPER(IC) 

CDP  =CDEP (IC) 

EDP-EDEP (IC) 

E  X  P  =  1  .  S 
BE  T  =  1  0/EXP 
KIN=CHF.ZY*SL  P**0 , 5 
QE-QP 

AE-  (  QE/K  IN  )  **EiET 

IF  ( AE  .  L.T  1  0E-5)  GO  TO  213 

VLL-QE/AF 

PPTH=CDP*AE**F ..DP 


I. Ml 


1 


WF.P=CPR*AE*#EPR 

RHR=AE/UEP 

HYR=DPTH 

UST“SQRT(32.2*HYR*SLP  ) 

IF ( IC . EQ . INLS )  GO  TO  241 
DELX=XSI(IC+1)-XSI(IC> 

GO  TO  242 

241  DELX=SLN-XSI < INLS) 

242  CONTINUE 
C 

C  COMPUTE  TRANSPORT  CAPACITIES  AND  TERM  DELTA  (EQ.  11) 

SUM  =  0 . 0 

DO  220  I  F  =  1 , NFR 
GL=GLAT ( I T , IF) 

CP=CI < IC , IF ) 

D  =  DMM (IF) 

DFT=DMM( IF ) /304 . 8 
U=VS( IF) 

IF <DMM( IF) , LE . 4  0 )  GO  TO  860 
IF  (  DMM(  IF  )  .  l.T  .  S  .  0  )GO  TO  861 
MEYERP=8 . 0 

CALL  MEYER (MEYERP , UST , RHB , DFT , GAMA , SNU , SPGR , WDTH , QE , VEL , CONC ) 
SCAP  < IF )=CONC*AE/l . 0E6/SPGR 
IF  <  DhM (IF)  GT.30.0)  SCAP<1F)=0.0 
GO  TO  862 

860  CALL  YANG ( DFT , UST , SNU ,  VEL , SLR , U , SPGR ,  CONC ) 

SCAP ( IF)=CCjNC*AE/1 . 0E6/SPCR 

GO  TO  862 

861  CALL  DUBOYS(D,DFT ,  AE , WDTH , SLP , SPGR , SNU , GAMA , CHI ,UST ,SFLOU> 
SCAP(IF) -S FLOW /VEL 

862  WATMAT  <  IF  )  ==CP*AE  +  GL*I)ELX/VEL 
IF (SCAP (IF)  EQ.0.0)  GO  TO  220 
5UM=SUM+WATMAT ( IF)/SCAP (IF) 

220  CONTINUE 
C 

C  COMPUTE  ACTIVE  (OR  ARMOR)  LAYER  THICKNESS 
P ARM=  0  0 
NARM=NFR 
DARM=DMM ( NFR ) 

DO  810  I F  =- 1  ,  NFR 
IF(SCAP(IF>  GT.0.0)  GO  TO  810 
IF( IF  EQ . 1 )  GO  TO  811 
IF(SCAPdF-i)  .GT.0.0)  CO  TO  811 
GO  TO  812 

811  NARM~ I F 
DARM=DMM ( IF) 

812  PARM=PARM  +  PC.E<C(  IC,  IF  ) 

810  CONTINUE 

IF ( FARM . CT  0.0)  GO  TO  814 

DO  81S  I FF  =  1 , NFR 

IF  A  =  NFR -I FF  f 1 

NARM=IFA 

DARM=DMM ( IFA) 

PARM  =  PCK<C( IC, IFA) 

IF ( PARM . CT . 0 . 0 >  GO  TO  814 
81S  CONTINUE 

814  AR MHT  = 10  0  (UDARM/PARM/304  8 

C  COMPUTE  VUL1JMF  OF  RED  MATERIAL  FRACTIONS  IN  ACTIVE  LAYER 
DO  813  IF =1, NFR 

RE DM A  1  ( IT  )  ARMHT*PCRC( I C , I F > *WDT H/ 1 0 0 . 0 


813 


IF<1. 0-SUM)  222,223,224 

continue; 

DEPOSITION  LOOP 
TBM-0  .  0 
BEDUP=0  0 
DO  230  J= 1 ,  NFR 
IF  =  NFR - J  + 1 
CP=CI < IC , IF ) 

IF(SCAP(IF)  .  EQ  .  0  0)  CO  TO  231 
RESCAP=SCAP < IF )*< 1 . 0-SUM ) 

DF.PO=0  0 

IF  CRESCAP  GE  0  0  )  GO  TO  232 
IFCABSCRESCAP  )  .  GT  UATMAI  (  IF)  )  C,0  TO  231 
DEPO=ABS(RESCAP ) 

GO  TO  232 
DEPO-WATMAT ( IF  ) 

BT-2  0*VS< 1F)*DTS/RHB 
IF (BT . L  T  1  0 )DFPO-BT *DEPO 
UJAT  MAT  (  IF  )  -  UATMAT  <  IF) -DEPP 
ADI)  =  0  .  0 

IF(SCAP(IF)  GT. 0  0)  ADD-  DEPO/SCAP ( IF ) 
SUM-SUM-ADD 

UPDATE  LOAD  AT  NODE  POINTS 
CAP -SCAP  < IF)/AE 
CCC=UATMAT< IF>/AE 
XKP=CAP-0  999*CP 
XK  C  =  C AP - 0  999*CXC 
ADD-0  .  0 

IF  (CAP  GT  0.0)  ADD'CLL  *f.AP  /  (  AE*XKP  *XKC  ) 

IF < ADD  LT  0.0)  ADD- 0  0 
FNLR=1  0+ADD 
DELT=DELX*FNLR/VEL 
IF  <  DELT . LT . DTS  >  GO  TO  220 U 
IFUC.EQ.  INLS)  GO  TO  233 
CCI-CI < IC+1 , IF) 

CC<IC+1,IF)  =  CCI+<CCC-CC1  )*i>TS/BFL  I 
GO  TO  234 

SEDIMENT  OUT  FLOW  FROM  CHANNEL  BLOCE: 

IF( IT  EQ  1)  CO  TO  23G 
CCI-G< IT-1 , 1 , IF )/0< IT-1  , 1  ) 

GO  TO  236 
CCI-CI < IC , IF) 

G  <  I  T  ,  I  ,  IF  )  -'  ( F.G I  +  (  CCC-CCI  )  *I)TS/I>f  L  T  )  #0  <  IT  ,  I  > 

GO  TO  234 

0  XDEL- DLLT * VKI../I  N!  R 

IF  < IC  EQ. 1 )  GO  TO  2201 
CC,l!  =  CC(IC,  IF) 

GO  TO  2202 

1  CCB~CUP ( IT , IF ) 

2  CC( IC+1 , IF)-CCB+ (CCC-CCEO*DEL X/XDEL 
IFUC.F.Q  I NLS  >  G<  IT  ,  I  ,  IF  >-CCt  JL  +  )  ,  IF  >*Q<  IT  ,  I  > 
BFDMAT ( IF ) -BEDMAT ( IF  >  +DEPO*DTS/Dr LT 

TBM  =  TBM  +  BE  DMA1 < IF  ) 

BEDUP  -BF.DUP  +  DEP O'# DTS  / DELT/UDTH/P OR  <  IF  ) 

CONTINUE 

CHANGE  IN  BED  ELEVATION  AND  BED  LAYER  COMPOSITION 
BEDEL V(  IC  )  =BE  DI.LV  <  IC)+BE  DUP 
DO  239  IF= 1 , NF  R 

PCBC  < IC, IF ) =BEDMAT  < I F ) / E PM* 1 0 0  .  0 


GO  TO  204 

223  CONTINUE 
C 

C  EQUILIBRIUM  LOOP 

DO  240  IF=  1 , NFR 
CP=CI  ( IC  ,  IF' ) 

CCC=WATNAT(IF)/AE 
CAP  =SCAP (IF>/AE 
XKP=CAP-0  999*CP 
XKC=CAP-0  9?9*CCC 
ADD=0  0 

IF (CAP  . GT . 0  .  0)  ADD=CEL*CAP/(AE#XKP#XKC> 

IF < ADD . LT  0  0)  ADD-0  0 
FNl.R  =  1  0+AI)D 
DELT+DELX*FNLR/VEL 
C  UPDATE  l OAD  AT  NODE  POINTS 

IF ( DELT  L  T  DTS )  CO  TO  2300 
IF<  IC  EQ  I NLS )  GO  TU  243 
CCI=C1 (IC+1 , IF) 

CC(  IC+1 , I F  > ~CC I  +  ( CCC-CCI )*DTS/DCLT 
GO  TO  240 

C  SEDIMENT  OUTFLOW  FROM  CHANNEL  BLOCK 
243  IF (IT  EQ  1 )  GO  TO  244 

CCI=G( IT-i , I , IF)/Q( IT-i  ,  I  ) 

GO  TO  24S 

2*4  CC I -C I ( I C , I F ) 

24S  G ( IT , I , IF )=(CC1+ (CCC-CCI) *DTS/DEL T ) *Q C IT ,1) 

GO  TO  240 

2300  XDEL- DLLT  fcUEL/FNLR 

IF ( IC  EQ. 1 )  GO  TO  2301 
CCB=CC ( I C , I F ) 

GO  TO  2302 

2301  CCB-CUP ( I T , I F ) 

2302  CC< IC+1 , IF)=CCB+(CCC-CCB>*DELX/XPEL 
IFUC.FQ  I  NLS)  G(IT,I,IF)=CC(IC  +  1,IF)#Q(IT,I) 

240  CONTINUE 
GO  TO  204 

224  CONTINUE 
C 

C  EROSION  AND  ARMORING  CALCULATIONS 
C 

C  COMPUTE  VOLUME  ENTRAINMENT  MATRIX 
DO  201  I  FA- 1  ,  NFR 
DEN=0  0 

LOOP-NFR-IFA+1 
DO  282  J  =  1 .LOOP 
IF-J+IFA-1 

282  DEN-DFN+COEF ( J ) *BEDMAT ( IF ) 

DO  283  J=1 .LOOP 
IF-J+IFA-1 

IF (DEN  EQ  00)  GO  TO  030 

ACCM ( IF , TFA)=COEF( J)*BEDMAT< IF  > /DEN*BEDMAT ( IF A ) 
GO  TO  283 

830  ACCM< IF , IFA)=0 . 0 

283  CONTINUE 
281  CONTINUE 
C 

C  EROSION  LOOP 

DO  270  IF  =1 , NFR 
TCA( IF ) =0 . 0 
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DO  279  IFA=i,IF 

279  TCA  ( IF)-TCA(IF)+ACCM(IF,IFA> 

278  CONTINUE 

DO  237  I  F  = 1 , MFR 
237  ERS< IF)=0 . 0 

DO  284  I F- 1 , NFR 

IF  <  SCAP (IF)  ,  EQ  0.0)  GO  TO  284 

RESCAP=SCAP ( IF >*( 1  0-SUH) 

I F  <  RESC AP  LE  O  0>  GO  TO  285 

IF  C  RESCdP  .GE.TCA(IF) )  GO  TO  288 

DO  831  J  =  1  ,  IF 

DO  287  Jl=--i  ,  IF 

IF  <  J1 . EQ . 1 )  I F  A  = I F 

IF(Ji.GT.l)  IFA-Ji-i 

IF  <  SC  AP  <  I  F  A  )  FQ  0  0)  GO  TO  287 

resc:ap=scap (  ifa)*(  i  o-sum> 

CNTBN= ACCM  (  IF  ,  IF  A) /FLOAT  (  IF  > 

IF  <  RESCAP  GT.CNIDN)  GO  TO  280 
ER  OSN-ERO*R  I.SCAP 
GO  TO  289 

288  EROSN^I  RO*CNT BN 

289  ERG  (  IFA)'ERS(  II  A)+EROF>N 
Sim=-SUM+EROSN/SCAP  (  IFA) 

IF (SUM  CF  1.0)  GO  10  288 

287  CONTINUE 
831  CONTINUE 

288  CONTINUE 

DO  290  IFA-  1  , 1 F 

IF (SCAP (IF  A)  IQ  0  0)  CO  TO  290 
EROSN  ERO*ACCM< IF ,  IFA  > 

ERS< IFA)=EKS( II  A )+FROON 
SUM-SUM t E R0SN/5CAP ( IFA) 

290  CONIINUI 

284  CONTIMUt 

C  UPDATE  LOAD  A  I  NODI  POINTS 

285  CONTINUE 
TBM  -  0  0 
DEDUN-O  0 

DO  2S0  IF-) .NFR 
CP  Cl  (  IC.  ,  I  F  ) 

UA1  MAI  (IF;  WA1MAT  (  IF  )  t[.RS(  IF  > 

CCC=UA1  MAT  (  I  F  )  /,*,£ 

CAP  =  SC  AH ( )t  )/AE 
XKP  -  C AP  -  0  9 9 9 ♦ C P 
XKC-CAP  0  999HCCC 
ADD- 0  0 

I F"  (  CAP  GT  .  0  0)  AI>P~CF  I  *CAF>/(AF.*/KP*XKC) 

IF  (ADD  L.  1  0  0  )  ADD-  0  0 

FNL.R  =  1  0t ADD 

DELT-DEl X#F  NLR/UEL 

IF  (DELT.L.T  DTS)  GO  TO  24  0  0 

IF(  IC  EQ  INLS)  GO  10  2S1 

CCI  =CI ( IG  +  1  ,  IF ) 

GC  (  IC  +  1  ,  IF  )  -CC I  +  <  CCC— CC I  >*DT3/I)FLT 
GO  TO  252 

C  SEDIMENT  OUTFLOW  FROM  CHANNEL  BLOCK 
251  IF (IT  EQ.l)  CO  TO  253 

CCI~G(IT-1,I,IF)/Q(IT-1,I) 

GO  TO  254 

253  CCI  =  CI  < IC , IF) 
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254  GUT,  I  , IF)=(CCI+<CCC-CCI>*DTS/DELT)*Q(IT,1> 

GO  TO  252 

2400  XDEL=DELT*VEL/FNLR 
IFUC  .  EQ . 1 >  GO  TO  2401 
CCB=CC(IC,IF) 

GO  TO  2402 

2401  CCB-CUP ( IT , IF ) 

2402  CCUC  +  l  ,IF)=CCB+(CCC-CCB)*DELX/XDEL 
IFUC.EQ.  INLS)  G  (  IT  ,  I ,  IF )  =CC  <  IC+ 1 ,  IF  >  *Q <  IT  ,  I  > 

252  TERS-ERS(IF)*DTS/DELT 

IF ( TERO . GT . BEDMAT ( IF) )  TERO=BEDMAT ( IF ) 

BE DM AT  ( IF ) “BEDMAT (IF ) -TERO 
TBM=TPM+BEDMAT( IF) 

C  UPDATE  BED  ELEVATION 

BEDWN-BEDWN+TERO/WDTH/PGR (IF) 

250  CONT I NUF 

BEDEL  V(  IC  )  BLDFI.V  (  IC  )  -BEDWN 
C  ADJUSTMENT  OF  ACTIVE  l.AYLR  THICKNESS 
IF ( TBM  GT  0  0)  CO  TO  294 
DO  293  IF  -  1  ,  NFR 

293  PCBC  < IC  ,  1 F ) =PCBI ( IC , IF  > 

GO  TO  204 

294  CONT 1 NUF 
P ARM-0  0 

DO  295  I F - 5  ,  NFR 

PCBC< IC , IF ) = BEDMAT  < IF ) /TBM* 1 00  0 
I F  <  SC.AP  (IF)  GT  0  0)  GO  TO  295 
PARM  =  PARM-t  PCBC(  IC  ,  IF  ) 

295  CONTINUE 

I F ( P AR M  Cl  0  0)  GO  TO  820 

DO  821  1  F F  = 1  .NFR 

IFA=NFR-IFF+1 

NARM  =  II  A 

DAR  M  =  DMM ( If  A) 

P ARM- PCBC ( IC . I F  A ) 

IF ( PARM  GT  0  0 )  GO  TO  820 
821  CONTINUE 

820  CUR ARM* 1  0  0  0 /? ARM*DAR M/30  4  8 

ADHT - CURARM- TBM/WDT  H 
IE ( ADHT  IE  0  0 )  GO  TO  204 
ADARM -ADHT 
TBM  -0  0 

DO  296  IF  -  1  .NFR 

BE. DMA T  (  IE  )  -  BFDMAT  (  IE  )  t  ADARM*PCB  I  (  IC  ,IF)/1  00  0*UDTH 

296  TBM=TBM« BEDMAT ( IF) 

DO  297  IF  -1,  NFR 

297  PCBC ( IC , IF ) - BEDMAT < I F ) /TBM* 10  0  .  0 
GO  TO  204 

213  CONTINUE 

C  DEPOSITION  CAUSED  BY  TERMINATION  OF  FLOW 


DO  260  IF  - 1  , NFR 
IFUC.EQ  I NLS )  GO  TO 

261 

261 

CCUC+l,  IF)  =  0.0 

GO  TO  262 

GUT, I,  I F )  =- 0  .  0 

262 

BEDMAT (IF) “BEDMAT ( IF 

)  +  CIUC,IF>*AE 

260 

CONTINUE 

204 

CONTINUE 

RETURN 

END 
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S  UBROUT INE  DUBO  YS  ( D  ,  DFT  ,  AE  ,  WDTH  ,  SLF  ,  SPGR  ,  SNU  ,  GAMA  ,  CHI , 

&UST,SFLOUI) 

THIS  SUBROUTINE  COMPUTES  THE  TRANSPORT  RATE  OF  NONCOHESIVE 
SEDIMENTS  USING  A  DUBOYS  TYPE  FORMULA ( GRAF ,  1971) 

CALL  SHIELD < DFT ,SPCR , SNU , GAMA , UST , TC > 

CHI=CHI/(D**0  7S ) 

TAO=GAMA*AE/WDTH*SlP 
SFLOU=WDTH*CHI*TAO* ( TAO-TC ) 

I F  <  SFLOU . LE .0.0 )SFLOW=0 . 0 

RETURN 

ENI) 

SUBROUTINE  MEYER < ME YERP , USTAR , RHP , D , WEIGHT , V  T  SC ,S , UI DTH , Q , V , GB , XO ) 

THIS  SUBROUTINE  COMPUTES  THE  TRANSPORT  RATE  OF  NONCOHESIVE 
SEDIMENTS  USING  THE  BEDLOAD  FORMULA  DEVELOPED  BY  MEYER-PETER 
AND  MULLER ( 1948 ) 

FCTN<F,RO,REY)-(i  /SQRT(F)  ) -2  .  *  At  OG1 0  <  RO  >  - 1  1-1  +  2  .  *ALOG10  <  1  +9  35* 
*RO/<REY*SQRT<F)  )  ) 

G  =  32 . 17 

REY-4  0*RHB*V/VISC 

R0=2 . *RHB/D 

R STAR  =  2  #D*USTAR/VISC 

IF  <  RSTAR  GT  70.0)  GO  TO  7 


F  1  =  0  OOfe 
F2=0  09 
1  =  0 

FCT1 =FCTN(F1 , RO , RE Y  ) 

F3=0  S*<F1+F2> 

FCT3=FCTN(F'3,R0,REY.) 

IF(FCTl*r  CT3) 1  ,5,2 
F2-+F3 
GO  TO  3 
F 1  =F  3 

IF ( APS ( F  2  -F  1  )  -  0  0001)5,5,4 
1  =  1+1 

IF  <  I  L  T  100)  GO  TO  F. 

GO  TO  7 
F  -F  3 
GO  TO  O 

F  =  <  1  0/(2  *  Al  OC 1  0  <  R  () )  +  1  14))**?. 

AK=V*SQRT (F/H  )/USTAP 
AK=0  75 

Y  =  ( UST AR*USTAR  )/(  (S-l  0  )*C.*D) 

Gi=MEYERP*WIDTH* (UST  AR**3  0 ) * ( WEIGHT /G)*(S/(S-1 . 0 ) > 
CALL  SHIELD  CD  ,S  ,  VISE.  ,UCIGHT  ,USTAR  ,TC> 

C  TSTAR=TC/< <S-1 . 0 )*WEIGHT*D> 

TST  AR  =  0 .047 
G2=AK**1 . 5- ( T  ST AR /Y ) 

IF  <  G2  .  LT  .  0  .  0  )  C.0  TO  9 
G?=G2**1 .5 
GB=G1*G2 

XO=  <  GB* 1 . 0E6)/<«*UEICHT) 

GO  TO  10 
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9 


10 

C 


GB  =  0  0 
XO=0  0 
RETURN 
END 


SUBROUTINE  YANGTDFT , UST , SNU , V , SLP , W , S , CONC > 


THIS  SUBROUTINE  COMPUTES  THE  TRANSPORT  RATE  OF  NONCOHESIVE  SEDIMENTS 
USING  THE  TOTAL  LOAD  FORMULA  DEVELOPED  BY  YANG (1973) 


D  =  DFT 

A=UST*D/SNU 

IF( A  GL  70.0)  GO  TO  7 

VCW=2  5/ <  AL OGi 0  <  A > -0 . 06  >  +  0 .66 

GO  TO  8 

VCW=2 . OS 

ESP=V*SLP/W~VCU*SLP 
IFIESP )  9,9,10 
CONC=0  0 
GO  TO  11 

0  Fi=S. 435-0  286*ALOGi0<U*D/SNU)-0  4S74AL0G 1 0 ( UST/U ) 

FP.-  l  .  799-0  409*ALOG1  0  <W*D/SNU>-0  31 4*ALDC1 0 ( UST/W ) 

F3~AL0G1 0 (ESP ) 

E=F1+F2*F3 

c=io  o**r 

CONC=C 

1  CONTINUE 

RETURN 
END 

SUBROUTINE  SHIELD <D ,S,VISC,WI ICHT ,USTAR ,TC) 

THIS  SUBROUTINE  COMPUTES  THE  CRITICAL  BED  SHEAR  STRESS  DERIVED 
FROM  SHIELDS'  FUNCTION 
R E Y  -  UST  AK*l)/VISC 
IF  <  RF Y  CT  10  0)  CO  TO  1 
TC--=0  0S*(S-1  .  0 >*WE1CHT *D/REY**0  4 
GO  TO  3 

IFIRtY  GT  500  (1  )  GO  TO  2 

TC*0  022* ( S- 1  0 ) *WEIGHT*D*RFY**0  16 

GO  TO  3 

TC---0  0  6  *  (  S  - 1  0)*UFICHT*D 

CONTINUE 

RE IURN 

END 

SUBROUTINE  SETUF.I.  (  D  ,U  ) 

THIS  SUBROUTINE  COMPUTES  THE  SETTl  1NG  VEIOCITY  OF  SEDIMENT  PARTICLES 
DIMENSION  A ( 2 , 1 1 ) 

DATA  A(1,1),A<J ,?),A(l,3),A(l,4),A(l,5),A(l,6>,A<l,?>,A<i,8>, 

4A< 1 , 9 ) , A< 1,10) ,A(1 , 11 )/ 

40.04,0  06,0.10,0  J0,0  40,0  80,1  50, 200, 300, 700, 10.  0(1/ 

DATA  A(2,1),A(2,2>,A<2,3>,A<2,4>,A<2,5>,A<2,6>,A(2,7>,A<2,8), 

4A  <  2 , 9 ) , A  <  2 , 1 0 ) , A  (  2 , 1 1 > / 

40.14,0  32,0 .76,2  20  ,5 .30 ,10  50 , 16.90 ,20 .30 ,25 .60 ,39.50 ,44 .00/ 

IF<D  LE  0  04)  GO  TO  20 
IF  <  D  GE . 10.0)  GO  TO  21 
GO  TO  22 

20  W-0 . 14/0 . 04*D 
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21 


GO  TO  100 

W=4 . 5/3. 0# ( D-i 0 . 0 ) +44 . 0 
GO  TO  100 
22  CONTINUE 

DO  10  1=1 ,11 

IF  <  A< 1 , 1 )  GT.D)  GO  TO  11 
Di  =  A< 1,1) 

U1=A<2,I ) 

GO  TO  10 

11  D2=A<1,I> 

U2=A (2,1) 

GO  TO  12 

10  CONTINUE 

12  W=<U2-Uli)*<D-Di)/<D2-Di)+Wi 

100  RETURN 

END 
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